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Abstract 

Low  thmst  propulsion  systems  such  as  electrodynamic  tethers  offer  a  fuel-efficient  means  to 
maneuver  satellites  to  new  orbits,  however  they  can  only  perform  such  maneuvers  when  they  are 
continuously  operated  for  a  long  time.  Such  long-term  maneuvers  occur  over  many  orbits  often  rendering 
short  time  scale  trajectory  optimization  methods  ineffective.  An  approach  to  multi-revolution,  long  time 
scale  optimal  control  of  an  electrodynamic  tether  is  investigated  for  a  tethered  satellite  system  in  Low  Earth 
Orbit  with  atmospheric  drag.  Control  is  assumed  to  be  periodic  over  several  orbits  since  under  the 
assumptions  of  a  nearly  circular  orbit,  periodic  control  yields  the  only  solution  that  significantly  contributes 
to  secular  changes  in  the  orbital  parameters.  The  optimal  control  problem  is  constructed  in  such  a  way  as 
to  maneuver  the  satellite  to  a  new  orbit  while  minimizing  a  cost  function  subject  to  the  constraints  of  the 
time-averaged  equations  of  motion  by  controlling  current  in  the  tether.  To  accurately  capture  the  tether 
orbital  dynamics,  libration  is  modeled  and  controlled  over  long  time  scales  in  a  similar  manner  to  the 
orbital  states.  Libration  is  addressed  in  two  parts;  equilibrium  and  stability  analysis,  and  control.  Libration 
equations  of  motion  are  derived  and  analyzed  to  provide  equilibrium  and  stability  criteria  that  define  the 
constraints  of  the  design.  A  new  libration  mean  square  state  is  introduced  and  constrained  to  maintain 
libration  within  an  acceptable  envelope  throughout  a  given  maneuver.  A  multiple  time  scale  approach  is 
used  to  capture  the  effects  of  the  Earth’s  rotating  tilted  magnetic  field.  Optimal  control  solutions  are 
achieved  using  a  pseudospectral  method  to  maneuver  an  electrodynamic  tether  to  new  orbits  over  long  time 
scales  while  managing  librational  motion  using  only  the  current  in  the  tether  wire. 


v 


To  my  wife  for  her  steadfast  support 


vi 


Acknowledgments 


I  thank  Mike  Ross  for  his  assistance  in  efficiently  transforming  math  problems  into  numerical  ones  that 
are  well-conditioned  for  numerical  optimization  and  Joe  Carroll  for  providing  insights  into  tether  design 
and  controls  that  helped  form  the  basis  for  using  long  time  scales  in  optimal  control  problems.  I  also 
appreciate  the  guidance  of  my  research  advisor  and  committee  chair  Bill  Wiesel.  I  am  especially  indebted 
to  Bill  Baker  for  the  many  hours  he  spent  helping  me  apply  the  mathematical  method  of  averaging  to 
optimal  control  problems. 


Robert  E  Stevens 


Table  of  Contents 


Abstract . v 

Table  of  Contents . viii 

List  of  Figures . x 

List  of  Symbols . xii 

I.  Introduction . 1 

II.  Literature  Review . 5 

III.  Optimal  Orbital  Maneuvering . 8 

Dynamic  Model .  10 

Constraints .  14 

Three  Optimal  Control  Problems  and  Their  Solutions .  16 

IV.  Multiple  Time  Scales  -  Modeling  Earth’s  Tilted  Magnetic  Dipole . 29 

Solution  to  an  Optimal  Control  Problem  Using  Multiple  Time  Scales . 35 

Summary . 38 

V.  Tether  Libration . 40 

Equilibrium  and  Stability . 40 

Demonstration  of  Attitude  Control  Using  Feedback  Linearization . 55 

Libration  Control  over  a  Long  Time  Scale . 60 

VI.  Summary,  Conclusions  and  Future  Work . 74 

Appendix  A:  Derivation  of  Libration  Equations  of  Motion . 77 

Kinetic  Energy . 78 

Potential  Energy . i 

The  Lagrangian  Equations  of  Motion . 84 


viii 


Non-Conservative  Generalized  Forces 


87 


Appendix  B:  Reduced  Mass  Derivation . 101 

Appendix  C:  Tether  Moment  Of  Inertia . 104 

Appendix  D:  Operational  Limitations . 105 

Preventing  Reentiy . 105 

Validity  of  a  Straight  Tether  Model . 109 

Appendix  E:  Tether  Tension  Curves . 114 

Appendix  F:  Scaling . 123 

Scaling  the  Time  Variable  for  Derivation  of  the  Averaged  State  Equations  of  Motion . 123 

Scaling  the  OCP  for  Well-Conditioned  Numerical  Solutions . 124 

Appendix  G:  Derivation  of  Averaged  Orbital  Element  Equations  of  Motion . 126 

Appendix  H:  Propagation  of  Libration . 129 

Appendix  I:  Reference  Synopses . 132 

References . 137 


IX 


List  of  Figures 


Figure  1.  Electrodynamic  Tether  Force  Model . 1 

Figure  2.  Optimal  Control  in  Fourier  Space . 9 

Figure  3.  Control  Solution  for  Maximum  Altitude  Maneuver  Using  32  Nodes . 19 

Figure  4.  Flamiltonian  Profile  for  Maximum  Altitude  with  Drag  Solution . 19 

Figure  5.  Maximum  Altitude  Maneuver  Trajectories.  Stars  indicate  DIDO  solution,  lines  indicate 

instantaneous  state  propagation  using  the  optimal  control . 20 

Figure  6.  Maximum  Inclination  Control  Solution  with  No  Drag . 21 

Figure  7.  Maximum  Inclination  Maneuver  Trajectory  without  Drag . 22 

Figure  8.  Maximum  Inclination  Control  Solution  with  Drag . 23 

Figure  9.  Maximum  Inclination  Maneuver  Trajectory  with  Drag . 24 

Figure  10.  Minimum  Time  Orbit  Change  Control  Solution  without  Drag . 26 

Figure  11.  Minimum  Time  Orbit  Change  Trajectory  without  Drag . 26 

Figure  12.  Minimum  Time  Orbit  Change  Control  Solution  with  Drag . 27 

Figure  13.  Minimum  Time  Orbit  Change  Trajectory  with  Drag . 27 

Figure  14.  Current  Control  in  Time  Clock  Time  Domain  for  the  First  16  Revolutions  for  Minimum  Time 

Maneuver  with  Drag . 28 

Figure  15.  Earth’s  Tilted  Magnetic  Dipole  Geometry . 29 

Figure  16.  Control  Profile  Using  Tilted  Dipole  Model . 36 

Figure  17.  State  Trajectories  Using  Tilted  Dipole  Model  with  Drag . 37 

Figure  18.  Optimal  Control  with  Drag,  but  No  Earth  Dipole  Tilt . 38 

Figure  19.  The  controller  model  breaks  down  when  Earth  magnetic  dipole  tilt  is  excluded  for  a  long  term 
orbit  transfer.  Stars  indicate  the  model-derived  altitude  trajectory;  line  indicates  propagated  altitude 

trajectory  in  a  rotating  tilted  dipole  magnetic  field . 39 

Figure  20.  Standard  MSIS  Atmosphere . 42 


x 


Figure  21.  EDT  Attitude  Geometry  defining  the  in-plane  and  out-of-plane  libration  angles  9  and  (j) , 

respectively . 43 

Figure  22.  In-plane  equilibrium  points  for  180  km  and  200  km  circular  orbits . 51 

Figure  23.  In-plane  equilibrium  points  for  250  km  and  300  km  circular  orbits . 52 

Figure  24.  Tether  Tension  Diagram . 53 

Figure  25.  EDT  Attitude  Control  Using  Feedback  Linearization . 59 

Figure  26.  Libration  Squared  Function  and  Envelope . 66 

Figure  27.  Control  for  Minimum  Time  Orbit  Change  with  Libration  Control,  No  Drag . 68 

Figure  28.  Minimum  Time  Orbit  Change  State  Trajectory,  No  Drag.  Stars  indicate  DIDO  derived  libration 

envelope;  lines  indicate  propagated  instantaneous  state . 69 

Figure  29.  Control  Profile  for  a  Minimum  Time  Orbit  Change  with  No  Libration  Control,  No  Drag . 70 

Figure  30.  Minimum  Time  Orbit  Change  State  Trajectory  without  Libration  Control . 71 

Figure  31.  Control  Profile  for  a  Minimum  Time  Orbit  Change  with  Libration  Control  and  Drag . 72 

Figure  32.  Minimum  Time  Orbit  Change  State  Trajectory  with  Drag . 73 

Figure  33.  Rotating  Frame  Coordinates . 77 

Figure  34.  Straight  Tether  Integration . 79 

Figure  35.  Position  of  Endmass  1 . 82 

Figure  36.  Body  Frame  and  Rotating  Frame . 85 

Figure  37.  Tether  Subject  to  Atmospheric  Drag . 88 

Figure  38.  Tether  Element  and  Drag  Geometry . 91 

Figure  39.  Electrodynamic  Tether  with  Lorenz  Force  Loading . 98 

Figure  40.  Tether  Mass  Distribution  Geometry . 101 

Figure  41.  Maximum  Lorenz  Force  to  Drag  Force  Ratio  at  the  Magnetic  Equator . 108 

Figure  42.  Tether  Curvature  Due  To  Lorenz  and  Aerodynamic  Force  Distribution . 110 

Figure  43.  Curved  Tether  Geometry . Ill 

Figure  44.  Matlab  ode23  Solution  to  Homogeneous  Equation  vs.  Exact  Solution . 129 


xi 


List  of  Symbols 


A 

4 

B 

B * 

cd 

D 

F« 

I 

4 

j 

L,  L 
L 

MnM 

Q, 

K 

T 

V 

a 

d, 

e 

f 

h 

ti 

i 

k 

mi 

n 

P 

r  ,r 

t 

u 

V,. 

v a 

X 

4 

e 

<t> 

r 

rm 

Pe 

P* 

P,n 


=  Jacobian  of  f 

=  area  of  component  i 

=  local  magnetic  flux  density  vector 

=  ballistic  coefficient 

=  coefficient  of  drag 

=  drag  force  vector 

=  force  vector  on  component  i 

=  instantaneous  control  current  in  electrodynamic  tether 
=  maximum  allowable  rms  control  current  in  tether 
=  specific  moment  of  inertia 
=  straight  tether  length  vector  and  magnitude 
=  Lagrangian  function 

=  effective  mass  of  end-body  i  and  total  system  mass 
=  generalized  force  i 
=  radius  of  the  Earth 

=  kinetic  energy  or  slow  time  scale  variable  (context  dependent) 

=  potential  energy 
=  semimajor  axis  of  the  orbit 
=  tether  diameter 
=  orbit  eccentricity 
=  state  vector  time  derivative 

=  altitude  of  system  center  of  mass  (COM);  vertical  rectangular  eccentricity  vector  coord 
=  atmospheric  scale  height 
=  orbit  inclination 

=  horizontal  rectangular  eccentricity  vector  coordinate  in  the  perifocal  frame 
=  mass  of  component  i 
=  orbital  mean  motion 
=  orbital  period 

=  position  vector  and  magnitude  from  Earth  center  to  system  COM 
=  time 

=  control  vector 
=  velocity  of  component  i 
=  velocity  of  air  with  respect  to  system  COM 
=  state  vector 

=  tether  linear  density  [mass/length] 

=  ratio  of  maximum  Lorenz  torque  to  gravity  gradient  torque;  scaling  parameter 
=  out-of-plane  libration  angle 
=  flight  path  angle 
=  Earth  magnetic  dipole  moment 
=  effective  reduced  system  mass 
=  Earth-satellite  gravitational  reduced  mass 
=  system  reduced  mass 
=  orbit  true  anomaly 


xii 


v 


6  =  in-plane  libration  angle 

p  =  atmospheric  density 

p,  Pj  =  straight  tether  unit  position  vector  from  ml  to  m2 ,  position  of  component  i  w.r.t.  COM 

t  =  tether  tension 

oo  =  argument  of  perigee 

0)(, ,  0Jn  =  angular  velocity  in  rotating  frame  and  inertial  frame  respectively 

Q  =  right  ascension  of  the  ascending  node 

v|/  =  basis  in  Fourier  space 

(//  =  tether  curve  angle 

( ’ )  =  time  derivative  unless  otherwise  specified 

( ' )  =  time  derivative  with  respect  to  rotating  frame 


xiii 


I.  Introduction 


With  increasing  dependence  of  government  missions,  scientific  exploration,  and  commercial  ventures 
on  spaceborne  payloads,  it  is  critical  to  have  the  right  satellite  over  the  right  place  at  the  right  time. 
Currently,  most  satellites  are  confined  to  Keplarian  orbits  that  reside  above  the  “reasonable”  atmosphere. 
Conventional  rockets  do  not  permit  a  satellite  to  orbit  at  lower  altitudes  where  atmospheric  drag  is  non- 
negligible  nor  do  they  usually  allow  large  orbit  adjustments  over  a  long  lifetime  since  these  scenarios 
would  require  a  prohibitive  amount  of  propellant.  However,  a  low  thrust  propulsion  system  requiring  little 
or  no  propellant  could  permit  station-keeping  at  lower  altitudes  and  even  provide  some  limited  orbital 
maneuvering  capabilities.  Electrodynamic  tethers  (EDTs)  in  low  Earth  orbit  offer  an  attractive  alternative 
to  conventional  satellites  that  use  propellant-based  propulsion  systems  because  the  thrusting  forces  are 
derived  using  the  Earth’s  geomagnetic  potential.  Electrodynamic  tethers  are  electrically  conductive  wires 
extending  between  two  or  more  subsatellites  and  when  a  current  is  passed  through  the  wire  in  the  Earth’s 
magnetic  field,  a  Lorenz  force  is  generated  perpendicular  to  both  the  current  direction  and  the  direction  of 
the  local  Earth  magnetic  field  lines.  A  two-ball  EDT,  defined  as  two  subsatellites  joined  by  a  conducting 
tether,  is  depicted  in  Figure  1  showing  how  the  Lorenz  force  generated  by  running  current  through  the  wire 
may  be  used  to  overcome  drag  and  maneuver  the  satellite  pair. 


o 

o 

o 
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Figure  1.  Eiectrodynamic  Tether  Force  Model 


The  force  magnitude  depends  on  the  current  /  ,  length  of  wire  and  the  wire  orientation  with  respect  to  the 
local  magnetic  field  according  to  Lorenz’s  law,  F  =  7L  x  B  ,  where  L  represents  the  length  vector  between 
the  satellite  pair  and  B  represents  the  local  Earth  magnetic  flux  density  vector.  Controlling  the  current  in 
the  wire  through  variable  resistance,  the  satellite  system  would  be  capable  of  maneuvering  to  new  orbits 
without  propellant,  albeit  at  a  slower  rate  than  traditional  maneuvering  rockets.  A  capability  such  as  this 
would  enable  space  missions  requiring  orbiting  sensors  at  extremely  low  altitudes  or  those  requiring 
frequent  repositioning  of  satellites  by  way  of  orbit  transfers. 

Because  of  the  low  thrust  provided  by  an  EDT  system,  an  orbit  transfer  requires  a  long  time  to  reach  a 
desired  orbit.  Obtaining  optimal  control  solutions  for  satellites  that  maneuver  for  a  long  time  can  be 
challenging  and  computationally  intensive  when  instantaneous  state  dynamics  and  controls  use  dynamics 
expressed  using  short  time  scales.  Williams  demonstrates  an  approach  to  optimal  control  using  non-linear 
perturbation  equations  of  motion  as  dynamic  constraints  and  solves  an  optimal  control  problem  by  direct 
transcription  using  Non-linear  Programming  (NLP)  software.  This  method  is  shown  to  be  effective  in 
determining  controls  that  execute  a  modest  orbital  maneuver  using  an  electrodynamic  tether  for  thrust; 
however  the  optimization  solver  required  hundreds  of  collocation  node  points  to  capture  all  the  small  state 
variations  that  occur  for  a  maneuver  that  only  takes  a  single  day.  Many  nodes  were  required  to  fully  depict 
the  instantaneous  states  and  control  that  exhibit  periodic  behavior.  Hundreds  of  collocation  nodes 
correspond  to  thousands  of  optimization  variables  and  constraints  for  the  NLP  solver  to  compute.  The 
number  of  nodes  and  computation  time  required  to  perform  the  optimization  over  long  periods  of  time  can 
be  difficult  or  impossible  to  achieve  using  the  short  time  scale  model  and  are  highly  susceptible  to  round¬ 
off  errors.  In  many  low-thrust  maneuvering  situations  the  instantaneous  orbit  state  will  vary  only  slightly 
from  Keplarian  motion  within  an  orbital  period  due  to  small  perturbations,  but  the  variations  tend  to  be 
periodic  in  the  short  term  and  cancel  out  over  the  long  term  leaving  only  slow  secular  state  changes. 
Addressing  long  term  behavior,  Carroll*  and  Tragesser  and  San  present  a  technique  of  non-optimal  periodic 
tether  control  that  uses  the  method  of  averaging  derived  from  perturbation  theory  enabling  control  of  the 
average  states  thus  avoiding  the  computational  burden  associated  with  controlling  the  rapidly  changing 
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instantaneous  states.  They  demonstrate  that  this  approach  is  good  for  determining  control  for  longer  time 
periods,  however  the  results  are  not  optimal  and  the  periodic  control  is  considered  to  be  unchanging 
throughout  the  trajectory.  Furthermore,  determining  the  control  requires  constraining  the  maximum  current 
which  is  less  straightforward  using  the  method  of  averaging  than  constraining  instantaneous  current  using 
the  short  time  scale  model.  The  problem,  therefore,  is  that  it  is  difficult  or  impossible  to  determine  optimal 
controls  for  an  EDT  performing  a  long  term  orbit  transfer  using  the  methods  of  control  currently  presented 
in  the  literature.  The  aim  of  this  dissertation  is  to  take  advantage  of  both  control  methods  to  achieve 
optimal  control  of  an  electrodynamic  tether  over  long  time  periods.  We  seek  to  modify  the  optimal  control 
problem  of  a  low  thrust  orbit  transfer  considering  that  we  already  know  something  about  the  dynamics  of 
the  system,  namely  that  it  is  nearly  periodic  when  the  EDT  is  continuously  thrusting.  Bearing  this  in  mind, 
we  may  dispose  of  the  dynamic  model  describing  the  rapidly  changing  instantaneous  behavior  in  favor  of  a 
dynamic  model  that  only  describes  the  secular  behavior  of  the  average  state  over  large  time  scales. 

The  research  objective  is  to  maneuver  an  electrodynamic  tether  to  a  new  orbit  over  many  revolutions  by 
posing  an  optimal  control  problem  in  the  context  of  large  time  scales  since  we  are  mainly  interested  in  the 
secular  behavior  and  not  the  periodic  behavior  occurring  during  each  revolution.  Although  this  research 
focuses  on  optimal  control  of  electrodynamic  tethers,  this  approach  to  optimal  controls  over  large  time 
scales  could,  in  principle,  apply  to  any  continuous  low-thrust  system. 

In  the  relevant  EDT  orbit  control  literature,  the  dynamic  models  used  are  limited  to  non-atmospheric 
environments  over  short  periods  of  time 1  or  they  ignore  attitude  dynamics  (libration)  and  long  term  optimal 
control2.  In  order  to  develop  a  real  system  that  will  operate  in  a  low  Earth  orbit,  drag  effects  and  libration 
must  be  addressed  and  either  included  in  the  controller  model  or  justifiably  ignored.  Because  the  source  of 
thrust  is  derived  from  the  Earth’s  tilted  magnetic  field,  it  is  also  important  to  include  the  effects  of  the 
Earth’s  rotation  on  the  satellite  motion.  To  achieve  the  primary  objective  of  determining  optimal  EDT  orbit 
transfers  spanning  many  orbits,  the  following  outline  describes  the  approach  taken  in  this  dissertation. 

•  Determine  optimal  long  term  maneuvers  for  nadir-pointing  tethers  ignoring  Earth’s  tilted 
magnetic  dipole  (Chapter  III) 

o  Develop  set  of  suitable  dynamic  equations  and  path  constraints 
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o  Pose  and  solve  optimal  control  problems,  first  ignoring  drag  and  then  including  drag 
o  Validate  optimal  control  solutions  by  propagating  with  a  “truth”  model 

•  introduce  the  effects  of  the  Earth’s  tilted  magnetic  dipole  into  the  dynamic  model  (Chapter  IV) 

o  Modify  the  dynamic  model  and  path  constraints  to  include  the  effects  of  a  tilted  Earth 
magnetic  field  that  rotates  once  per  day 
o  Pose  and  solve  optimal  control  problems  with  and  without  drag 
o  Validate  optimal  control  solutions  by  propagating  with  a  “truth”  model 

•  Introduce  tether  libration  into  the  dynamic  model  (Chapter  V) 

o  Perform  stability  analysis  of  tether  libration 

o  Modify  the  dynamic  model  and  path  constraints  to  include  the  effects  of  tether 
libration 

o  Pose  and  solve  optimal  control  problems  with  and  without  drag 
o  Validate  optimal  control  solutions  by  propagating  with  a  “truth”  model 
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II.  Literature  Review 


The  main  body  of  spacecraft  tether  literature  as  it  relates  to  this  research  may  be  divided  into  several 
categories.  In  the  area  of  tether  dynamic  analysis  and  control,  the  chief  motions  studied  are  orbital, 
librational  and  vibrational.  Models  used  for  these  motions  depend  on  the  application,  but  will  sometimes 
include  either  electrodynamic  forces  [Refs  1-2,  6-15]  or  aerodynamic  drag  [Refs  3,4  and  5],  but  very  few 
include  both.  Researchers  have  investigated  tether  control  strategies  using  tether  length  variation,  end-body 
drag,  thrusters,  and  in  the  case  of  electrodynamic  tethers  (EDTs),  wire  current  (references  provided  in 
forthcoming  discussion).  This  research  will  focus  on  the  unexplored  area  of  long  term  optimal  control  of 
an  EDT. 

The  first  category  of  relevant  literature  addresses  orbital  maneuvering  using  EDTs.  Most  of  these 
papers  discuss  system  design  issues  but  do  not  detail  controller  design.  Tragesser  and  San  describe  various 
EDT  current  controllers,  but  they  are  non-optimal  and  librational  motion  is  ignored.6  In  the  area  of  optimal 
control,  the  published  works  are  very  limited.  The  most  relevant  paper  discussing  the  optimal  control  of 
electrodynamic  thrusting  tethers  is  one  from  Williams. 7  The  dynamic  model  in  this  paper  ignores  the 
atmosphere  and  the  librations  are  not  explicitly  bounded,  however  the  paper  showcases  an  example  optimal 
orbital  maneuver  useful  for  comparison.  In  other  related  orbital  maneuvering  works,  long  term  EDT 
thrusting  strategies  were  published  first  by  Carroll*  and  then  Tragesser  and  San  for  no-drag  orbits,  however 
optimal  control  over  large  time  scales  is  not  addressed. 

There  is  a  large  body  of  work  addressing  the  second  category,  EDT  libration  analysis  and  control.  In 
the  case  of  electrodynamic  tether  models,  Pelaez  et.  al.  explore  the  stability  of  these  systems  assuming  a 
constant  tether  current  for  inclined 910  and  elliptical  orbits.  In  more  elaborate  analyses,  two  bar  tether 
models  were  employed. 1 1  Although  many  of  these  papers  do  not  address  control,  they  do  provide  insights 
into  the  behavior  of  unthrottled  active  electrodynamic  tethers  in  a  non-atmospheric  environment. 

Reference  12  shows  the  librational  instability  that  occurs  with  a  constant  current  EDT,  thus  control  is 
necessary  to  compensate  for  drag  while  simultaneously  maintaining  stability..  Without  control,  an  EDT 
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system  would  need  a  “self-balanced”  design  to  maintain  stable  attitude  dynamics  according  to  Ref  13 .  Ref 
14  concludes  that  EDT  control  can  be  employed  to  manage  instabilities  for  orbits  with  eccentricity  less 
than  0.35.  Hoyt  presents  a  method  to  stabilize  using  feedback  control. 15  There  are  other  methods  of 
attitude  control  for  hanging  EDTs  besides  using  torques  due  to  Lorenz  forces  about  the  center  of  mass 
(COM).  Williams  describes  a  method  of  libration  stability  control  using  tether  length  variation16,  as  does 
Yu  for  orbits  with  e<0.3,  however  only  in-plane  motion  is  considered.17  Thrusters  have  also  been  proposed 
for  libration  control, ls  however  by  using  propellant,  this  method  defeats  the  stated  purpose  of  using  EDTs 
in  the  first  place.  Most  controller  designs  used  one  of  the  linear  techniques  like  Ref  19  which  describes 
thruster  and  tension  control  using  LQR  methods.  Some  papers,  such  as  Refs  20  and  21 ,  present  nonlinear 
control  methods  (feedback  linearization)  to  maneuver  between  equilibrium  points.  A  combination  of 
control  methods  is  presented  in  Ref  22  where  both  electromagnetic  forces  and  length  rate  are  used  to 
manage  librations. 

De  Matteis  and  De  Socio  caution  against  instabilities  due  to  atmospheric  density  gradients  in  very  long 
tethers  (>75  km)  that  could  lead  to  a  destabilizing  libration  resonance  at  altitudes  lower  than  240  km.2j 
The  culprit  in  this  case  was  that  a  long  librating  tether  would  be  subject  to  very  different  drag  forces 
throughout  large  pendular  swings.  However  for  the  tether  lengths,  operational  altitudes  and  allowable 
librations  considered  in  this  work,  the  density  variations  are  relatively  minor  and  this  effect  is  ignored. 

Another  category  of  tether  literature  is  devoted  to  vibration  and  mode  shapes.  Von  Flotow  shows  that 
a  tether  under  the  uniform  loading  of  an  electrodynamic  and  aerodynamic  force  will  tend  to  sag  in  the 
middle  with  a  slow  first  lateral  mode  of  vibration  (slow  relative  to  the  longitudinal  vibration).  Using  the 
fact  that  the  period  of  the  first  lateral  mode  of  vibration  is  long,  he  approximates  the  tether  to  be  in  a  state 
of  quasi-static  equilibrium  in  the  shape  of  a  section  of  a  circular  arc.24  This  shape  and  vibration 
approximation  is  used  in  determining  the  maximum  current  limits  for  the  system  (Appendix  D).  Other 
vibration-related  work  includes  control  of  an  electrodynamic  tether  through  input  shaping  to  reduce 
vibrations  and  librations25  and  vibrations  due  to  a  constant  Lorenz  force  load.26  Watanabe  suggests  a  bang- 
bang  current  control  providing  input  shaping  to  reduce  vibrations  and  librations  while  thrusting  with  an 
EDT.27  Williams  investigates  control  of  flexible  tethers  using  electromagnetic  forces  and  a  movable 
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attachment.28  These  analyses,  however,  do  not  include  atmospheric  effects.  De  Matteis29  presents 
equations  of  motion  that  include  aerodynamic  effects  in  modeling  vibrational  behavior  of  non¬ 
electrodynamic  tethers.  There  are  several  authors  who  develop  controllers  for  non-electrodynamic  tethers 
using  tension  control  or  reels  for  length  variation. 30  Others  focus  on  the  deployment  and  retrieval  phases.  31 

The  remainder  of  the  tether  literature  is  largely  aimed  at  specific  design  studies  or  missions.  Many 
authors  have  addressed  designing  tethers  to  operate  efficiently  and  safely  in  the  space  environment.  Bare 
wires  efficiently  collect  electrons  to  produce  the  current  used  for  thrust.32,33  Porous  tapes  have  been 
proposed  and  investigated  to  increase  the  survivability  where  micrometeors  can  sever  a  thin  wire  and  end  a 
mission.34,35  Several  other  papers  were  written  to  support  SEDS  (Small  Expendable  Deployer  System)  and 
other  specific  space  demonstrations.36,37,38  Estes  et.  al.  document  lessons  learned  from  the  various 
missions  that  have  deployed  in  space.39  A  good  reference  covering  all  the  general  topics  reviewed  is  the 
book  Dynamics  of  Space  Tether  Systems  by  Beletsky  and  Levin.40 

In  the  area  of  optimal  control  of  low  thrust  multirevolution  transfers,  Ross,  Gong  and  Sekavat 
propose  a  technique  that  manages  the  high  frequency  content  of  optimal  solutions.  Solutions  are  achieved 
by  solving  a  large  time  scale  optimal  control  problem  using  a  small  number  of  nodes.  Applying  Bellman’s 
principle,  they  then  iteratively  solve  the  problem  and  propagate  the  control  solution  along  smaller  sections 
of  the  original  optimal  path,  thus  capturing  all  its  detailed  high  frequency  components.  This  general 
method  has  the  advantage  of  solving  large  time  scale  optimal  control  problems  while  still  avoiding  the 
aliasing  common  when  there  are  not  enough  collocation  node  points  to  resolve  the  high  frequency  content. 
A  more  exhaustive  list  of  the  literature  reviewed  along  with  a  brief  synopsis  of  their  pertinent  contents  is 
included  in  Appendix  I. 

The  electrodynamic  tether  literature  provides  ample  coverage  of  libration  stability  analysis  and  control, 
however  only  a  few  papers  address  orbital  maneuvering.  Williams  paper  on  optimal  orbit  transfer  and 
Tragesser  and  San’s  maneuvering  approach  using  periodic  control  stood  out  as  the  two  works  most  relevant 
to  the  research  presented  in  this  dissertation.  Starting  with  key  concepts  extracted  from  these  two  papers,  a 
new  approach  to  optimal  orbital  maneuvering  is  presented. 
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III.  Optimal  Orbital  Maneuvering 


In  this  chapter,  we  examine  optimal  maneuvers  using  a  two-ball  EDT  as  defined  in  the 
introduction  (Figure  1).  Because  of  the  low  altitudes  considered,  the  trajectories  account  for  atmospheric 
drag  and  are  nearly  circular  therefore  the  orbital  equations  of  motion  may  be  expanded  about  the  very  small 
eccentricity.  This  assumption  is  good  for  orbits  with  small  eccentricities  as  long  as  errors  remain  within  the 
tolerance  of  the  spectral  algorithm  used  for  optimization.41'42  Furthermore,  the  maneuvers  are  known  to 
occur  over  many  orbital  revolutions,  so  the  small  oscillatory  changes  in  the  orbital  parameters  that  are 
evident  over  short  time  scales  (within  each  revolution)  are  averaged  out  leaving  only  the  secular  changes 
that  occur  over  long  time  durations  (many  revolutions).  See  Appendix  F  for  a  full  discussion  on  different 
time  scales.  The  only  control  we  have  at  our  disposal  to  perform  the  desired  maneuvers  is  the  current  in  the 
wire  using  variable  resistance.  The  dynamic  behavior  of  the  EDT  in  its  slowly  changing  orbit  is 
predictably  periodic,  consisting  of  a  linear  combination  of  sinusoidal  functions  of  true  anomaly,  or 
equivalently,  time.  Because  motion  of  a  constantly  thrusting  EDT  deviates  from  Keplarian  motion  in  a 
periodic  manner  over  each  orbit,  a  controller  that  is  also  periodic  over  the  orbit  will  contribute  to  secular 
changes  in  the  orbital  parameters  over  a  long  time  as  shown  in  Appendix  G.  Other  contributions  of  the 
controller  are  averaged  out  in  the  long-term.  Therefore  we  assume  periodic  control  current,  /  ,  modeled 
using  the  relevant  terms  of  a  Fourier  series. 

/(t',u(r))  =  Im  (w,  ( T)  +  u2  (T)cosv  +  m3  (T)sinv  +u4  (T)cos2v  +u5  (7T)sin2v/)  (1) 

where  v  is  the  true  anomaly  and  Im  is  the  maximum  allowable  rms  control  current.  To  highlight  the  fact 
that  the  controlled  Fourier  coefficients  vary  only  over  large  time  scales,  we  write  them  as  functions  of  T . 
The  slow  time  scale  variable  T  is  a  scaled  version  of  the  clock  time,  t ,  which  is  itself  proportional  to  the 
true  anomaly.  All  state  variables  and  controls  that  are  functions  of  T  change  very  slowly  and  are 
considered  constant  over  short  time  intervals.  For  a  more  complete  discussion  on  time  scaling,  see  Chapter 
V  and  Appendix  F.  The  control  in  Fourier  space  with  bases  Y(V)  =  [l,cosv,sinv,cos2v,sin2v]r  is 
therefore  completely  defined  by 

u  (r)  =  [zq ,  u2 ,  m3  ,  u4  ,u5f 
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(2) 


such  that  the  current  in  the  clock  time  domain  is  given  in  Eq.  (1)  by  /  =  /mYr  (Vju  ( T ) .  With  the  control 
written  in  this  form,  the  approach  to  optimal  control  is  viewed  in  the  Fourier  space  where  the  goal  is  to 
determine  the  time  dependent  Fourier  coefficients,  u  ( T ) .  that  minimize  a  given  cost  function  for  a 
trajectory  subject  to  the  time -averaged  dynamic  equations  of  motion.  A  pseudospectral  method  of  dynamic 
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Figure  2.  Optimal  Control  in  Fourier  Space 


optimization  is  employed  using  DIDO  software41'42  to  solve  the  subject  optimal  control  problems  yielding 
the  optimal  control  coefficients  and  path  discretized  over  large  periods  of  time.  The  diagram  in  Figure  2 
shows  that  optimal  control  problems  exhibiting  periodic  behavior  may  be  transformed  into  a  Fourier  space 
using  the  method  of  averaging.  This  eliminates  dependence  on  fast  time  variables  ( v,  t )  and  the  resulting 
averaged  states  and  controls  will  only  depend  on  the  slow  time  variable  ( T ).  The  problem  posed  in  this 
Fourier  space  may  now  be  solved  by  an  optimizer  producing  average  states  and  controls  that  change  slowly 
over  time. 

If  it  is  desired  to  capture  the  instantaneous  states  for  subsections  of  the  averaged  trajectory,  then 
the  general  method  of  multirevolution  optimal  control  proposed  by  Gong,  Ross  and  Sekhavat  would  be 
suitable.  In  this  way,  instantaneous  optimal  controls  could  be  determined  for  sections  of  the  optimal 
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averaged  path.  This  method  was  tested  over  a  shorter  time  span,  where  an  optimal  control  problem  was 
posed  and  solved  using  this  antialiasing  algorithm  that  applies  Bellman’s  principle  to  capture  high 
frequency  content  of  an  optimal  trajectory.  The  algorithm  was  used  to  solve  a  fixed  time  (6  orbits), 
maximum  inclination  problem  (constrained  to  use  positive  current  only)  initially  using  32  nodes.  Then  the 
algorithm  iteratively  solved  subsections  of  the  optimal  path  using  32  nodes,  propagated  the  resulting 
optimal  control  which  provided  a  new  initial  condition  for  the  optimizer  until  the  desired  end  conditions 
were  reached.  The  final  solution  used  16  iterations  and  was  able  to  capture  fast  time  periodic  behavior  of 
some  of  the  states  {0,(/),a  ),  but  the  secular  behavior  was  difficult  to  observe  due  to  the  few  number  of 
revolutions.  Longer  term  maneuvers  were  not  attempted  due  to  the  number  of  iterations  required  to  capture 
all  the  fast  time  dynamics  of  the  instantaneous  states.  The  method  did  demonstrate,  however,  that  it  could 
potentially  be  used  to  solve  for  sections  of  the  optimal  averaged  path  without  placing  any  assumptions  on 
the  controller  (i.e.  not  required  to  be  periodic)  perhaps  taking  advantage  of  some  of  the  higher  order  effects 
availed  by  using  instantaneous  state  dynamics.  The  remainder  of  this  dissertation  will  focus  on  the 
construction  of  and  solution  to  optimal  EDT  control  problems  in  Fourier  space,  and  will  begin  by  deriving 
a  dynamic  model. 

Dynamic  Model 

Orbital  changes  due  to  the  relatively  weak  Lorenz  forces  generated  along  the  electrodynamic  tether 
occur  over  many  orbital  revolutions.  The  EDT  is  modeled  as  a  “dumbbell”  consisting  of  two  end  bodies 
tethered  together  with  a  taut  (i.e.  positive  tension)  4  km  copper  wire.  The  Lorenz  force  generated  along  the 
wire  containing  electric  current  is  given  by 

F  =  TLxB  (3) 

where  /  represents  tether  current  (the  control),  B  represents  the  Earth’s  local  magnetic  flux  density  vector, 
and  L  is  the  tether  length  vector  pointing  in  the  direction  from  the  upper  end  mass  to  the  lower  one  (see 
Appendix  A).  The  tether  geometry  and  current  direction  that  yields  a  positive  transverse  thrust  is  shown  in 
Figure  1. 

The  local  magnetic  flux  density  for  an  Earth-orbiting  satellite  is  modeled  as 
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(4) 


-2sin(<y  +  v/)sin/ 

cos  (Vo +  v)  sin/ 

= 

B, 

cos; 

e 

Bn 

where  ym  represents  the  Earth’s  magnetic  dipole  moment,  i  is  the  inclination  relative  to  the  magnetic 
equator  and  Br ,  Bt ,  and  Bn  represent  the  magnetic  flux  density  vector  components  in  the  radial,  transverse 
and  orbit  normal  directions  respectively  (i.e.  er ,  et  and  en  directions).  At  the  equator,  a  force  of  0.1  N 
distributed  along  a  one  amp,  4  km  EDT  is  achievable  at  an  altitude  of  270  km,  which  can  be  the  same  order 
of  magnitude  as  the  atmospheric  drag  at  that  altitude  depending  on  the  physical  characteristics  of  the  tether 
and  end  bodies.  To  ensure  the  satellite  orbits  longer  than  a  few  days,  the  control  system  will  need  to  apply 
a  constant  average  current  in  order  to  provide  constant  in-track  thrust  that  will  compensate  for  drag  forces 
acting  in  the  opposite  direction.  The  problem  of  drag  compensation  is  exacerbated  when  the  EDT  orbits  at 
a  higher  inclination  since  the  out-of-plane  component  of  the  magnetic  field  which  produces  the  required  in¬ 
track  thrust  is  reduced  (see  Eq.  (4)).  Drag  magnitude  depends  on  the  physical  properties  and  dimensions  of 
the  EDT,  the  atmospheric  density  and  satellite  velocity.  For  a  near  circular  orbit,  the  drag  force  on  the 
entire  tether  system  is  given  by 

D  =  (5) 

2  r 

where  p{r)  represents  the  average  air  density  at  radial  distance  r,and  B  is  the  average  ballistic 


C  A 

coefficient  of  the  entire  tether.  Flere  the  ballistic  coefficient  is  defined  as  B *  =  — - —  where  C.  is  the 


average  coefficient  of  drag,  A  is  the  average  cross-sectional  area  perpendicular  to  the  velocity  vector,  and 
m  is  the  system  mass.  Modeling  the  atmosphere  as  an  exponentially  decaying  density  using  a  scale  height 
h  ,  we  can  expand  about  the  small  eccentricity  and  approximate  the  average  density  through  first  order  as43 


P{r)  =  P0e  '•  ~p{a)  1  + 


aecosv 


h 


(6) 


where  the  radial  distance  has  been  approximated  as  r  «  a(l-ecosv) . 
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Gravity  gradient  torque  tends  to  keep  the  tether  nadir-pointing  with  libration  that  is  assumed  to  be  small 
(libration  is  addressed  in  Chapter  V),  so  the  acceleration  due  to  the  Lorenz  force  in  Eq.  (3)  is  given  by 


F  = - ^-(cos/et  -  cos  (v  +  <u)  sin /en )  (7) 

mr  v  ' 

Recognizing  that  the  orbits  of  interest  at  this  low  altitude  are  nearly  circular,  we  ignore  ()(  e2 )  and 

higher  order  terms  and  write  the  equations  of  variation  for  the  five  classical  orbital  elements  as 

da  2a,  r  „ 

_  F  +  D  -et 

dt  ai¬ 


de  1 


( 


dt  na~e 
dco  1 


-r  (F  +  D)-et 


1  +  - 


1 


dt  nae\  1  +  ecosv 
di  rcos(v'  +  ®) 


.  ,  „  rcot;sin(v  +  o) 

sin  v  (F  +  D)-e, - ' - -F  e„ 


na 


(8) 


dt  na ~ 

dCl  rsin(v  +  o) 


-F  -e„ 


F  •  e„ 


dt  «a2sin/ 

where  n  is  the  mean  motion  of  the  satellite  (Ref  43  pp.  84-85).  Expanding  these  equations  of  motion  about 
the  small  eccentricity  using  r  k  ~  a  k  (l  +  fecosv)  and  ignoring  second  and  higher  order  terms,  we  write 

the  general  perturbation  equations  of  motion  for  a  nadir  pointing  tether  in  terms  of  the  true  anomaly.  This 
is  the  only  variable  that  changes  significantly  on  a  short  time  scale. 


—  »  2Cacos/7(v)(l  +  4<?cosv')-2D  1  +  f  2  +  -^rlecosvl 

dt  \  \  h  J  ) 


de  _  . T/  , /  _  r  2  \  2D 

—  ~  C  cos// (v  2 cos v  +  5e cos’  v+e) - 

dt  K  A  ’a 


cosv  +  e  +  |  l  +  ^v  |ecos2v 
h 


/(v)(l  +  2e  cosv)  +  sini' 


C  cos  /  _  /  \  /  \  D 

- /(v)(2  +  5ecosv) - 

e  ae 


2  +  |  1  +  ——  |ecosv 

h 


dco  Ccos/sin(2o  +  2v) 
dt  2 

—  «  -Csin/7(v')cos2  (v  +  o)(l  +  2ecosv') 
dt 

dQ.  Cl(v) 

- « - -sin(2v  +  2o)(l  +  2ecosv) 

dt  2 

Ly  B'jup(a) 

We  have  let  C  = - ^  represent  the  term  resulting  from  thrust,  and  D  = - - — ^  represent  the  drag  term. 

nma  2  na 


(9) 


In  this  form,  these  equations  could  serve  as  dynamic  constraints  in  posing  our  optimal  control  problems, 
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however  due  to  the  rapid  variation  of  true  anomaly  with  each  revolution  we  would  need  to  discretize  the 
problem  with  enough  node  points  for  the  solver  to  capture  the  motion  of  each  varying  element  with  each 
revolution.  This  is  the  approach  Williams  used  (Ref.  1)  to  achieve  optimal  control  solutions  for  short  time 
scale  problems.  Since  we  are  only  interested  in  the  secular  state  changes  of  the  EDT  orbit  over  long  time 
scales,  we  use  the  method  of  averaging  to  eliminate  the  small  oscillations  that  occur  within  each  revolution 
which  effectively  approximates  the  nonautonomous  system  in  Eq.  (9)  as  an  autonomous  averaged  one.44 
This  is  achieved  by  recognizing  that 

dt  =  —  \l  -  2e  cosv  +  0[e2y^dv  (10) 

and  then  integrating  over  2ttN  ( N  =1,2,...).  Because  the  average  states  vary  slowly  with  time,  they  are 
considered  constant  over  the  short  time  periods  of  integration  and  are  removed  from  the  integrand.  The 
fast-time  variable,  v  ,  always  appears  in  the  argument  of  a  sine  or  a  cosine  function,  therefore  integrating 
Eq.  (9)  with  respect  to  v  will  yield  non-zero  results  only  when  the  control  current,  /  ,  is  itself  periodic  (i.e. 
it  is  a  combination  of  sine  and  cosine  functions  of  v  ).  A  current  that  is  purely  dc  will  produce  secular 
motion  in  semi-major  axis  and  inclination  since  the  first  two  derivatives  in  Eq.  (9)  would  yield  non-zero 
values  after  integration.  Because  an  EDT  depends  on  the  Earth’s  magnetic  field  for  propulsion,  the  orbits 
of  interest  remain  very  close  to  the  Earth  and  are  therefore  nearly  circular.  To  avoid  singularity  neare  =  0 , 
we  will  substitute  two  equinoctial  coordinates  for  the  eccentricity  and  argument  of  perigee  in  Eq.  (9).  The 
new  coordinates  are  the  eccentricity  vector  components  defined  as  h  =  esin  co  and  k  =  ecos  co  .  The 
average  state  equations  of  motion  are  derived  in  Appendix  G  using  the  periodic  current  defined  in  Eq.  (1), 
and  are  written  as 
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Secular  changes  to  the  orbit  state  are  now  expressed  over  a  large  time  scale,  AT  = - .  The  state  vector 

n 

x  now  represents  the  average  orbital  state  values  rather  than  the  instantaneous  values  and  is  written  using  a 
quasi  equinoctial  element  set,  i.e.  x (T)  =  \a,  h,  k,  i,  Q]r  .  Notice  that  these  average  states  vary  slowly  over 

long  time  scales  (indicated  by  T )  and  are  considered  constant  “within"  each  revolution.  The  average  state 
equations  of  motion  are  thus  devoid  of  the  short  time  scale  variable,  true  anomaly.  From  the  first  equation 
in  (1 1)  we  see  that  the  average  drag  effect  due  to  the  air  density  (in  the  drag  term,  D  )  primarily  affects  the 
average  change  in  semi-major  axis.  To  a  lesser  extent  drag  decreases  the  h  and  k  states  and  has  a 

circularizing  effect  since  e  =  yfh2  +k2  .  With  the  secular  equations  of  motion  in  hand,  we  now  turn  to 
constraining  the  allowable  tether  current  to  values  that  are  within  the  system  power  limitations. 


Constraints 

To  determine  the  optimal  controls  for  the  system  described  by  Eq.  (11),  we  need  to  solve  for  the 
periodic  control  coefficients,  u  (  T  ) .  Besides  enforcing  the  initial  state  conditions  as  event  constraints,  the 

control  current  must  also  be  bound  to  remain  within  an  available  power  limit  which  is  itself  defined  by  the 
electron  collection  capabilities,  ohmic  losses,  voltage  current  and  other  factors.  For  a  description  of 
electron  collection  in  the  ionosphere  and  the  associated  limitations  see  Ref  45 .  Because  the  control  in  Eq. 
(1)  is  defined  using  the  rapidly  changing  true  anomaly  we  cannot  simply  bound  the  instantaneous  periodic 
current  between  a  minimum  and  maximum  value  since  we  need  to  keep  our  averaged  equations  of  motion 
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devoid  of  short  time  scale  variables.  To  properly  bound  the  control  then,  we  need  to  define  a  path 
constraint  that  is  a  function  of  the  slowly  varying  Fourier  control  coefficients,  u  (/’ )  .  The  approach  used 
herein  limits  the  average  power  available  for  thrust  which  in  turn  places  bounds  the  on  the  rms  current.  For 
a  given  constant  wire  resistance  R  and  average  power  limit,  PavgMax  ,  the  maximum  allowable  rms  current  is 

defined  by  Joule’s  law  combined  with  Ohm’s  law 

(12) 

K 

The  actual  electric  current  rms  value  over  one  orbit  (period)  is  defined  by 

lL=^jl2(»,v)dv  (13) 

in  0 

For  the  periodic  current,  this  value  is  (using  Eq.  (1)) 

IL  =  II  (j',2  +  |(«2  +  «3  +  «4  +  «S  )]  (14) 

Using  Eq.  ( 14)  we  may  express  the  path  constraints  in  terms  of  the  controls.  The  path  constraint  for  the 
control  is  written  as 

sMT))=IL-P^fL^  (15) 

which  places  an  upper  bound  on  the  rms  control  current  throughout  the  transfer.  Choosing  a  proper  value 
for  the  maximum  allowable  rms  control  current  is  addressed  in  Appendix  D. 

This  path  constraint  approach  has  the  double  advantage  of  averaging  out  any  parameters  periodic  with 
the  orbit  that  affect  the  available  thrust  current,  such  as  diurnally  varying  ionospheric  electron  density,  as 
well  as  eliminating  the  short-time  variable,  the  true  anomaly.  The  event  constraints  (constraints  on  states 
at  specific  times  during  the  trajectory)  are  comprised  only  of  the  initial  conditions  and  are  written  as 

e(x(T0 ))  =  [a(T0 ),  h(T0  ),k(T0),i(T0jf  (1 6) 

Finally,  states,  controls  and  time  are  bounded  by  upper  and  lower  limits  (denoted  using  subscripts  ‘u’  and 
T  respectively).  These  box  constraints  are  written  as 
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(17) 


x,  <x(r)<x„ 
u,  <u(r)<u„ 

T  <T  <T 

10l  —  10  —  10u 

T  <T  <T 

±fl  -  ./  -  1fu 

Now  all  the  pieces  are  in  place  to  construct  and  solve  optimal  control  problems  that  will  maneuver  an  EDT 
to  a  new  orbit  over  many  revolutions  while  overcoming  drag  by  controlling  nothing  but  current  in  a  wire. 


Three  Optimal  Control  Problems  and  Their  Solutions 

Three  sample  maneuvers  were  chosen  to  demonstrate  large  time  scale  optimal  control  because  of  their  slow 
secular  orbital  changes  that  occur  over  many  revolutions.  The  tether  modeled  in  all  three  problems  is  4  km 
long  and  2  mm  in  diameter  (based  on  TiPs,  a  nonconducting  tether  system  deployed  in  1996).  The  system 
mass  and  average  cross-sectional  area  is  500  kg  and  8  m2  respectively.  The  first  two  problems  outline  the 
optimal  control  problem  setup,  solution  and  results  for  maximizing  the  average  altitude  and  inclination,  and 
serve  as  benchmark  problems  since  other  authors  have  investigated  similar  non-optimal  problems.1 2  The 
third  problem  provides  an  example  optimal  control  problem  and  solution  that  achieves  a  minimum  time 
orbit  change  occurring  over  500  revolutions  using  only  40  nodes  in  the  discretized  optimization  problem. 
All  problems  were  solved  using  DIDO,  an  optimization  software  package  that  discretizes  and  solves 
general  optimization  problems  using  a  pseudospectral  method.46  Even  though  the  derivation  that  produced 
Eq.  (11)  required  integration  over  a  hypothetical  integer  number  of  revolutions,  the  optimizer  does  not  need 
to  discretize  the  trajectory  over  the  same  integer  multiples  since  the  dynamic  equations  of  motion  are 
established  for  the  continuous  average  state,  not  the  instantaneous  state.  This  average  state,  however,  is 
meaningless  unless  the  total  maneuver  time  is  long  enough  to  span  several  periods.  This  akin  to  ensuring  a 
sample  interval  is  big  enough  to  capture  all  the  desired  frequencies  in  signal  analysis. 

Verification  of  the  optimal  control  solution  was  achieved  by  evaluating  the  Hamiltonian  output  by 
DIDO.  To  demonstrate  the  accuracy  of  the  model  used  as  the  dynamic  constraint  in  these  problems,  the 
output  Fourier  coefficient  controls  were  converted  into  the  time  domain  and  then  used  to  propagate 
instantaneous  states  using  Eq.  (8). 
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Maximum  Final  Altitude 


Consider  the  scenario  where  there  is  a  need  to  tow  an  object  (spacecraft,  debris,  etc.)  to  a  higher  orbit  in 
the  same  orbital  plane  using  an  EDT.  For  the  sake  of  testing  the  algorithm  against  a  known  solution  we 
seek  the  maximum  altitude  an  EDT  can  reach  in  50  orbital  revolutions  with  no  drag.  In  this  case  we  expect 
that  a  direct  current  in  the  nadir-pointing  tether  will  provide  maximum  thrust  in  the  direction  of  the  velocity 
to  spiral  the  spacecraft  out  to  a  higher  orbit. 2S  Although  we  may  actually  want  to  control  the  other  orbital 
elements  to  a  desired  end  state,  we  seek  only  this  known  solution  for  this  benchmark  problem.  The  optimal 
control  problem  is  written  as  the  following. 

Minimize  Cost:  J  =  -af 

Subject  to: 

Dynamic  Constraints  x  (T)  =  f  ( x  (T) ,  u  (7° )) 

Event  Constraints  e(x(T0))  =  [6648  km, 0, 0.001, 30°]7 

Path  Constraints  g,  (u  (T ))  =  I2rms  -  2.25  <  0  Amps2 

where  x(j)  is  the  average  state  change  and  f  (x(T),u(T))  =  Ax/A T  .  Box  constraints  in  Eq.  (17)  are 
also  enforced  where  we  have  chosen  the  bounds  to  be 

x„  =[16000  km,0.4,0.4,80\180°]r 

x,  =[6638  km, -0.4, -0.4,15°, -180“]r 

u„  =[l,V2,V2,V2,V2]r  (18) 

11/  =  -u 

T0=  0 

Tf  =  50  P 

where  P  is  the  orbital  period  at  T  =  0 .  The  initial  states  h  and  k  correspond  to  an  eccentricity  of 
0.00 1  and  an  initial  argument  of  perigee  of  zero.  Before  using  the  optimization  solver,  the  states  and  time 
were  scaled  to  span  values  of  order  1  to  make  the  problem  numerically  well-conditioned.41  47  Solving  the 
problem  using  DIDO  yields  the  control  history  shown  at  the  top  of  Figure  3,  and  the  bottom  of  the  figure 
shows  the  control  transformed  into  the  short  time  scale  domain,  in  this  case  just  a  direct  current.  The 
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average  altitude  and  inclination  trajectories  are  shown  in  Figure  5,  where  the  stars  indicate  the  DIDO 
solution  at  distinct  times  (spanning  large  time  scale  steps)  and  the  lines  indicate  the  propagation  of  the 
instantaneous  state  values  using  DIDO  derived  controls  and  a  Matlab '  ode  solver.  As  expected,  in  order  to 
perform  a  maximum  climb  maneuver  the  solution  indicates  that  the  controller  should  drive  a  maximum 
allowable  direct  current  through  the  wire  to  accomplish  the  large  transverse  thrust  needed  to  boost  the  orbit. 
Starting  at  an  altitude  of  270  km,  this  EDT  can  climb  about  130  km  in  about  3  days  without  drag. 
Introducing  drag  into  the  dynamic  constraints  does  not  affect  the  control  profile,  but  reduces  the  achievable 
altitude  change  in  the  given  number  of  orbit  periods  (50)  to  about  117  km.  In  reality  we  would  need  to 
contend  with  libration  control  and,  at  times,  adverse  battery  conditions  that  could  limit  power  available  for 
tether  thrusting.  However,  in  principle,  modest  maneuvers  can  be  accomplished  if  they  are  not  time 
critical.  Because  there  is  no  explicit  time  dependence  in  the  Lagrangian  of  the  Hamiltonian  of  this  optimal 
control  problem  (Eq.  (19)),  the  resulting  Hamiltonian  should  be  constant,  i.e.  H  =  0  .  The  Lagrangian  of 
the  Hamiltonian  is 

/7  =  //  +  iugg,+H^x+n>  (19) 

where  the  Hamiltonian  is  given  by  H  =  krf  and  /,  represents  the  costates.  The  covector  functions 
associated  with  the  path  constraint,  state-variable  box  constraints  and  control-variable  constraints  are 
represented  by  /ug  ,  p  and  p  respectively.  DIDO  uses  the  Covector  Mapping  Principle48  to  produce 

adjoints  and  the  Hamiltonian  as  part  of  the  solution.  To  check  optimality  the  output  Hamiltonian  was 
plotted  and  it  was  revealed  that  it  was  indeed  constant  as  shown  in  Figure  4. 
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Max  Climb  Fourier  Control  Coefficients,  u(T) 
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Figure  3.  Control  Solution  for  Maximum  Altitude  Maneuver  Using  32  Nodes 


Max  Climb  Solution  with  32  Nodes 


Figure  4.  Maximum  Altitude  Maneuver  Trajectories.  Stars  indicate  DIDO  solution;  lines 
indicate  instantaneous  state  propagation  using  the  optimal  control. 
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Figure  5.  Hamiltonian  Profile  for  Maximum  Altitude  with  Drag 


Maximum  Final  Inclination 

From  Eq.  (1 1)  it  is  evident  that  a  carefully  and  constantly  applied  dc  control  current  could  indeed 
compensate  precisely  for  drag  to  maintain  altitude,  however  it  would  come  at  the  expense  of  a  secular 
decay  of  the  inclination  after  a  long  time,  which  may  be  undesirable.  To  maximize  the  final  inclination 
achievable  in  a  fixed  time  (now  for  500  revolutions),  we  write  the  same  optimal  control  problem  as  in  the 
previous  example  with  the  following  exceptions. 

Minimize  cost:  J  =  —i. 

Subject  to: 

7}  =  500P 

e(x(r0))  =  [6648  km, 0, 0.01, 30°  f 
g2(x(r))  =  /r  +  /i2-eo2=  0 

where  the  new  path  constraint,  g2  (x(7’)j ,  ensures  a  constant  eccentricity  transfer. 
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As  a  test  case,  we  first  look  for  the  no  drag  solution  (i.e.  atmospheric  density  terms  in  Eq.  (11)  are  zero), 
then  compare  it  with  the  solution  that  accounts  for  drag.  The  32-node  DIDO  control  solution  to  the  no  drag 
problem  is  depicted  in  Figure  6. 

The  contrast  between  the  two  plots  of  the  same  control  in  Figure  6  clearly  shows  the  advantage  of 
solving  the  optimal  control  problem  using  Fourier  coefficients  over  a  large  time  scale.  Attempts  to 
discretize  and  optimize  this  control  problem  using  instantaneous  states  and  their  respective  dynamic 
equations  of  motion  (Eq.  (8))  for  this  long  term  trajectory  would  require  thousands  of  nodes  and  run  the 
risk  of  round-off  errors  and  long  solution  times.  Propagating  the  instantaneous  states  using  the  optimal 
control  output  produces  the  trajectory  shown  in  Figure  7  where  the  magnified  inserts  clearly  show  the 
instantaneous  fast  time  dynamic  behavior. 

Because  there  is  no  drag  to  contend  with,  the  optimal  solution  indicates  that  it  is  best  to  mainly  use  an 
ac  current  that  has  double  the  orbital  frequency,  i.e.  a  combination  of  m4  and  u5  within  constraints.  This 
result  is  consistent  with  Refs  2  and  8  which  indicate  that  to  achieve  maximum  inclination  change  the 
control  strategy  is  to  drive  a  current  such  that  /  =  —J21m  cos2(  v  +  .  Flere,  it  is  assumed  that  the  path 

constraint  in  Eq.  (15)  is  active  which  bounds  the  peak  amplitude  of  this  ac  input  to  V2 Im .  Transforming 
this  result  into  the  Fourier  coefficient  controller  described  in  Eq.  (1)  we  see  that  the  control  solution  is  the 
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Figure  6.  Maximum  Inclination  Control  Solution  with  No  Drag 
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same,  only  expressed  in  the  context  of  the  partial  equinoctial  set.  To  achieve  a  maximum  final  inclination, 


the  control  may  be  written 


I 


-V2/m 


cos(2v  +  2o)  =  ~\/2 4,  (cos 2® cos 2v  -  sin  2® sin  2v) 

( k2-h 2  „  2 kh  .  „  \ 

- 2 —  cos  2v - y~  sin  2v 


(20) 


In  this  form  we  recognize  the  Fourier  coefficients  for  the  second  mode  cosine  and  sine  functions  as 


=  -V2/ 


(  1,2 


k-h 


2  A 


(21) 


Max  Inclination,  No  Drag  Solution  Using  32  Nodes 


Figure  7.  Maximum  Inclination  Maneuver  Trajectory  without  Drag 
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The  optimal  controls  calculated  using  Eq.  (21)  are  consistent  with  the  corresponding  control  Fourier 
coefficients  determined  by  DIDO  (within  an  error  2-norm  of  0.04).  This  trajectory  uses  some  negative  dc 
thrust  to  decrease  altitude  while  increasing  the  orbit’s  inclination.  The  Fourier  control  coefficients 
displayed  in  Figure  6  show  that  the  tether  current  controller  initially  uses  a  small  negative  dc  component  to 
descend  to  the  lowest  allowable  altitude  in  order  to  maximize  the  final  inclination.  Controlling  the 
spacecraft  in  such  a  way  increases  the  orbit  inclination  from  30°  to  31.19°  in  about  a  month.  This  strategy 
outperforms  a  similar  constant  altitude  maneuver  by  0.04°.  When  drag  is  considered,  the  control  strategy  is 
altogether  different  because  more  of  the  limited  available  current  must  be  constant  dc  in  order  to 
compensate  for  the  increased  drag  as  seen  by  comparing  Figure  6  and  Figure  8.  We  see  from  Eq.  (1 1)  that 
a  large  positive  dc  coefficient  tends  to  reduce  the  inclination.  There  is  a  penalty  for  orbiting  where  the 
atmospheric  density  is  higher  because  more  power  is  expended  simply  to  maintain  altitude  which  causes 
inclination  to  decay  and  less  power  available  to  maximize  the  inclination.  In  this  case,  the  strategy  is  to 
climb  to  a  lower  density  altitude,  level  off  to  increase  inclination  then  descend  again  to  the  minimum 
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Figure  8.  Maximum  Inclination  Control  Solution  with  Drag 

allowable  altitude  taking  advantage  of  the  largest  possible  inclination  gain  opportunities  as  shown  in  Figure 
9.  The  initial  climb  comes  at  the  expense  of  inclination  gain,  however  overall  the  satellite  achieves 
maximum  inclination  change  because  it  operates  in  a  lower  average  drag  environment  and  does  not  need  to 
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Max  Inclination  Solution  with  Drag  Using  32  Nodes 
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Figure  9.  Maximum  Inclination  Maneuver  Trajectory  with  Drag 

expend  as  much  power  to  maintain  altitude.  After  a  month  of  thrusting  in  a  reduced  average  drag 
environment,  the  satellite  achieves  an  inclination  gain  of  1 .0°  outperforming  a  constant  altitude  maneuver 
by  0.25°.  Because  this  maneuver  occurs  over  so  many  revolutions,  it  would  be  near  impossible  for  short 
time  scale  optimization  to  yield  a  solution  to  this  problem.  The  problem  is  complex  when  attempting  to 
solve  in  the  clock  time  domain  but  it  is  reduced  to  a  simple  Zermelo  problem*  in  Fourier  space.  The  next 
example  problem  will  demonstrate  how  to  apply  this  method  to  a  more  general  orbit  transfer. 


In  1923  German  mathematician  Ernst  Zennelo  posed  the  problem  of  navigating  a  boat  from  point  A  to 
point  B  in  minimum  time  factoring  in  wind  and  current.  The  solution  is  not  a  straight  line  path.  Ref.  Jean 
van  Heijenoort,  1967.  From  Frege  to  Godel:  A  Source  Book  in  Mathematical  Logic,  1879-1931.  Harvard 
Univ.  Press. 
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Minimum  Time  Orbit  Change 


Having  looked  at  the  baseline  orbital  maneuvers,  we  now  turn  our  attention  to  determining  the  controls 
for  a  minimum  time  orbit  change  involving  a  desired  final  altitude  and  inclination  while  maintaining  a 
constant  eccentricity.  In  this  example  we  start  by  using  our  initial  states  from  the  first  example  and  then 
construct  the  optimal  control  problem  to  achieve  a  10  km  climb  and  a  one  degree  inclination  increase, 
while  maintaining  a  constant  eccentricity  of  0.005,  in  the  quickest  time.  We  write  the  problem  as 

Minimize  Cost:  J  =tf 

Subject  to: 

x(r)  =  f(x(r),u(r)) 
e0  (x(r0))  =  [6648  km, 0, 0.005, 30°  f 
ef(x(7}))  =  [a/,t/]r  =[6658  km,31°f 
ft  (u(r))  =  /L-  2.25  <0  Amps2 
g2(x(T))  =  h2+k2-e20=  0 

Box  constraints  are  still  those  listed  in  relations  (17)  and  Eqs.  (18),  but  since  this  problem  has  a  free 
final  time,  we  write 

T0+5<Tf  <5xl04H 

The  control  solution  without  drag,  depicted  in  Figure  10,  indicates  that  the  strategy  is  to  initially  apply 
a  negative  dc  control  current,  indicated  by  w, ,  to  descend.  The  controller  needs  to  apply  large  ac  control 
components  cycling  at  twice  the  orbital  frequency  to  reach  the  desired  inclination  (i.e.  large 
u4  and  u5  components),  all  while  avoiding  large  components  cycling  at  the  orbital  frequency,  namely 
u2  and  u3 ,  which  are  large  contributors  to  eccentricity  change.  The  dc  component  is  nearly  zero  for  the 
majority  of  the  trajectory  and  then  reverts  to  positive  flow  at  the  end  of  the  trajectory  to  climb  to  the  final 
desired  orbital  altitude  (Figure  11).  When  drag  is  considered,  the  dc  component  of  the  control  current  is 
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throttled  (see  Figure  12)  such  that  the  satellite  initially  climbs  and  then  descends  to  the  final  orbit 
minimizing  the  cost  due  to  increased  drag  at  lower  altitudes  as  much  as  possible  as  shown  in  Figure  13. 
Contending  with  drag,  this  EDT  takes  an  additional  four  days  to  complete  the  maneuver. 


Minimum  Time  Orbit  Change,  Control  Fourier  Coefficients,  u(T),  No  Drag 
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Figure  10.  Minimum  Time  Orbit  Change  Control  Solution  without  Drag 


Minimum  Time  Solution  with  40  Nodes,  No  Drag 
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Figure  11.  Minimum  Time  Orbit  Change  Trajectory  without  Drag 


Controlling  the  slowly  varying  current  Fourier  coefficients  over  many  revolutions  has  the  advantage  of 
solving  long-term  problems  with  relatively  few  nodes  in  the  optimization  algorithm.  A  similar  problem 
solved  using  a  small  time  scale  and  exact  equations  of  motion  would  yield  the  instantaneous  states  during 
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each  revolution,  however  it  would  require  an  exorbitant  number  of  nodes  over  the  same  time  frame  to 
arrive  at  a  meaningful  solution.  The  periodic  current  would  require  at  least  four  nodes  per  orbit  revolution 
in  the  short  time  scale  domain  to  establish  a  control  current  that  avoids  the  node  points  aliasing  undesired 
harmonics.  The  first  day  alone  in  this  example  consists  of  32  control  current  cycles  (Figure  14)  which 
would  require  at  least  64  nodes  to  adequately  capture  all  the  cycles.  Using  large  time  scales  and  averaged 
states,  we  have  solved  a  multirevolution  orbital  maneuvering  problem  using  40  optimization  nodes 
contrasted  to  the  two  thousand  nodes  that  would  have  been  required  using  a  short  time  scale  and 
instantaneous  states.  These  examples  demonstrate  that  solving  optimal  control  problems  in  Fourier  space 


Control  Fourier  Coefficients,  u(T),  Minimum  Time  Maneuver  with  Drag 

2, - T - , - , - 1 - 1 - 1 - T - ! - 

C/3 


Revolutions 


Control  in  Time  Domain 


.4I _ 1 _ 1 _ 1 _ i _ 1 _ 1 _ i _ i _ 1 _ 

0  50  100  150  200  250  300  350  400  450  500 

Revolutions 


Figure  12.  Minimum  Time  Orbit  Change  Control  Solution  with  Drag 


Minimum  Time  Maneuver  Solution  with  Drag  using  40  Nodes 
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Figure  13.  Minimum  Time  Orbit  Change  Trajectory  with  Drag 
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using  large  time  scales  and  time-averaged  orbital  states  has  significant  advantages  when  the  desire  is  to 
control  the  secular  behavior  of  a  continuously  operating,  low  thrust  satellite  system  over  a  long  time  rather 
than  the  instantaneous  periodic  behavior.  In  satellite  control,  a  rapidly  changing  periodic  variable  may  be 
averaged  out  leaving  only  the  dynamics  of  the  slowly  changing  variables.  In  this  dissertation,  a  method  of 
constructing  and  solving  a  large  time  scale  optimal  control  problem  using  an  electrodynamic  tether  to 
maneuver  to  a  desired  orbit  has  been  investigated.  Optimal  controls  for  three  sample  maneuvers  spanning 
up  to  500  orbital  revolutions  were  determined  using  30  to  40  optimization  nodes  instead  of  the  hundreds  or 
thousands  of  nodes  required  using  the  instantaneous  clock  time  dynamics.  The  remainder  of  this 
dissertation  will  use  the  concepts  introduced  here  to  improve  the  controller  dynamic  model  by  including  the 
effects  of  the  Earth’s  rotating  tilted  magnetic  dipole  (Chapter  IV)  and  tether  libration  (Chapter  V). 
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Figure  14.  Current  Control  in  Time  Clock  Time  Domain  for  the  First  16  Revolutions  for 
Minimum  Time  Maneuver  with  Drag 
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Table  1.  Summary  of  Results 


Maneuver  Type 

No  Drag 

With  Drag 

Maximum  Final  Altitude 

Aa  =  1 30  km 

Aa  =  1 17  km 

Maximum  Final  Inclination 

Ai  =  1.19° 

© 

II 

< 

Minimum  Time  Orbit  Transfer 

Av  =  432  revs 

Av  =  499  revs 
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IV.  Multiple  Time  Scales  -  Modeling  Earth’s  Tilted 


Magnetic  Dipole 

Because  electrodynamic  tethers  depend  on  the  Earth’s  magnetic  field  to  generate  a  thrusting  force,  an 
accurate  model  of  this  field  is  required  to  accurately  control  the  spacecraft.  The  models  used  so  far  have 
assumed  that  the  Earth’s  magnetic  dipole  moment  vector  is  aligned  with  the  Earth’s  poles  and  the  magnetic 
field  remains  constant  as  the  EDT  orbits  through  it.  In  reality  the  Earth’s  magnetic  dipole  moment  vector  is 
tilted  with  respect  to  the  North  Pole  by  about  1 1.5  degrees  (according  to  NASA)  and  rotates  with  the  Earth 
once  per  day.  Since  the  local  magnetic  field  vector  at  any  given  point  in  the  EDT’s  orbit  cycles  with  a 
period  of  one  sidereal  day,  the  controller  must  account  for  this  effect  in  the  model.  Fortunately  this  effect  is 
predictably  periodic  and  may  be  included  in  our  existing  model  of  averaged  state  dynamics.  The  diagram 
in  Figure  15  depicts  the  planes  containing  the  geographic  equator,  magnetic  equator,  and  the  EDT  orbit 
where  the  magnetic  equatorial  plane  rotates  about  the  North  geographic  pole  vector  (N). 

The  inclination  and  argument  of  latitude  at  epoch  of  the  satellite  with  respect  to  the  magnetic  equatorial 
plane  are  represented  by  im  and  am  respectively  for  a  dipole  that  is  tilted  by  5  .  The  argument  of  latitude  at 


N 


Figure  15.  Earth’s  Tilted  Magnetic  Dipole  Geometry 
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epoch  of  the  satellite  with  respect  to  the  geographic  and  equatorial  planes  are  given  respectively  by 

a  =  co  +  v 
°cm=a-y 

The  inclinations  with  respect  to  the  two  equatorial  reference  planes  are  related  using  the  law  of  cosines  of 
spherical  trigonometry. 

cos  im  =  cos  8  cos  i  +  sin  8  sin  i  cos  (Q  -  Q . ) 

S'1'  im  =  ^~cos2im 

where  we  assume  the  satellite  to  be  in  an  orbit  such  that  0  <  im  <  90°  Vt  >  0  and  Q  is  the  angle  from  an 
inertial  reference  direction  to  the  intersection  of  the  two  equatorial  planes  in  the  direction  of  the  longitude 
of  the  ascending  node  in  the  geographic  equatorial  plane.  This  angle  varies  with  time  over  a  medium  time 
scale  (i.e.  a  day)  and  is  related  to  the  true  anomaly  by 

>1 

where  77  is  a  scaling  factor.  Applying  the  spherical  trigonometric  laws  of  sines  and  cosines  once  again,  we 

obtain  the  relationships  between  the  arguments  of  latitude  in  both  reference  equatorial  planes. 

sin  y  sin  im  =  sin  8  sin  (Q  -  Cle ) 
cos  y  sin  i  sin  im  =  (cos  8  -  cos  i  cos  im ) 

The  argument  of  latitude  at  epoch  in  the  magnetic  equatorial  plane  is  then  written 

sin  am  =  sin  («  -  / )  =  sin  a  cos  y  -  cos  a  sin  y 
cos  am  =  cos  (or  -y)  =  cos  a  cos  y  +  sin  or  sin  y 

The  local  magnetic  field  vector  direction  is  a  function  of  the  inclination  and  argument  of  latitude  at  epoch 
with  respect  to  the  magnetic  equator  may  now  be  expressed  in  terms  of  a  time  varying  function  that 
depends  on  the  same  orbital  parameters  referenced  to  the  geographic  equator.  In  the  satellite  frame  (Figure 
33  in  Appendix  A),  the  magnetic  field  vector  is 
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-2 


B  =  % 
r 


-2  sin  am  sin  im 
cos  am  sin  im 
cos  i 


sin  a 
sin/ 
cos  a 


sin; 


: — 7  (cos  8  -  cos  i  cos  im )  -  cos  a  sin  <5  sin  (Q  -  Qe ) 
(cos  8  -  cos  i cos  im )  +  sin  a  sin  8 sin  (Q  -  fle ) 


cos;,,. 


After  making  the  appropriate  substitutions  and  grouping  terms,  a  form  for  the  magnetic  field  in  terms  of  the 
orbital  elements  referenced  to  the  geographic  equatorial  plane  is  derived  as 

B  =  Bj  cosct+B^  sin  8  (22) 


where 


B, 


-2  sin  a  sin ; 
cos  a  sin ;' 
cos; 


B„ 


2  [sin  a  cos ;'  cos  (Q  -  Qc )  +  cos  a  sin  (Q  -  )] 

-  cos  a  cos ;'  cos  (Q  -  Qf )  +  sin  a  sin  (Q  -  Q.e ) 
sin;'cos(f2-fle) 


(23) 


The  dynamic  model  of  the  local  magnetic  flux  density  vector  may  be  decoupled  into  two  terms,  each 
affecting  the  motion  of  the  EDT  over  a  different  time  scale.  The  first  magnetic  field  term  B,  is  derived 
from  the  same  non-tilted  dipole  moment  model  used  in  Eq.  (4)  and  it  is  periodic  on  a  short  time  scale  of 
one  orbit,  but  does  not  change  with  respect  to  the  medium  time  scale  (i.e.  a  day).  It  may  be  averaged  over  a 
single  period.  The  second  magnetic  field  term  Bn ,  however,  is  periodic  on  a  medium  time  scale  of  // 

orbits  (//  >  1)  since  it  contains  terms.  For  optimal  control  problems  spanning  times  such  that 
v  »  2 nr]  this  magnetic  field  term  may  be  averaged  over  i]  orbits.  As  expected,  when  the  model  assumes 
that  there  is  no  tilted  dipole  moment,  then  Eq.  (22)  reduces  to  the  standard  model  used  in  Eq.  (4)  that  is 
periodic  over  one  period. 

The  Lorenz  force  due  to  a  control  current  driven  through  the  tether  is  given  as 
F  =  /LxB  =  /Lx  (Bj  cos  <5  +  B^  sin  (S')  =  F,  cos  +  Fn  sin;? 
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where  F,  represents  the  contribution  to  the  electrodynamic  force  that  is  periodic  over  a  single  orbit  and  Fn 

represents  the  contribution  that  is  periodic  only  over  77  orbits,  i.e. 

Fj  =7LxB, 

F„  =/LxBn 

The  F,  contribution  is  the  same  as  that  derived  in  Eq.  (7)  and  may  be  averaged  over  one  orbit  as  shown  in 
the  example  problems  in  the  previous  section  using  the  non-tilted  dipole  model.  The  Fn  term  contributes 

to  the  total  Lorenz  force  when  there  is  a  non-zero  tilt  in  the  Earth’s  magnetic  dipole  moment  and  the  vector 
direction  cycles  with  period  iTir/ .  For  a  nadir-pointing  EDT  using  the  satellite  frame  defined  in  Figure  33, 
this  force  contribution  is 

F,  =  a((B,-e.)et-(B11.et)e.) 

with  components 


F„er=0 

F^  -et  =  Cl  sin/(cos£TcosQc  +  sinQsinQf)  (24) 

Fn  en  =  Cl  [cos  a  cos;  (cos  Q  cos  Q(,  +sinQsinQe)-sina(sinQcosQ(,  -cos£TsinQe)] 

Ly 

where  C  = - ^  .  Assuming  that  states  and  controls  change  very  slowly  over  77  orbits,  the  average  states 

7777777 


change  due  to  this  force  contribution  only  when  there  is  a  control  current  that  is  resonant  with  Q  =  — .  We 

77 

will  therefore  define  a  control  current  as  in  Eq.(l)  that  now  includes  Fourier  control  coefficients 
m6  and  it-j  that  correspond  with  this  harmonic  as 

I  (v,u(T),?i(T))  =  7  m,  +  772  cos v  +  77 3  sin  V  +  77 4  cos2v  +  775  sin2v  +  776  cos  — +  77,  sin—  (25) 

\  U  r/) 

Substituting  Eqs.  (24)  and  (25)  into  the  perturbation  equations  of  motion  (Eq.(8))  and  changing  the 


independent  variable  to  the  true  anomaly  using  the  approximation  dt  *  —  (l  -2ecos  v)dv  we  write  the 


averaged  perturbation  equation  for  the  semi-major  axis  as 

Aa  =  Aa,  cos  5  +  A  a  sin  <7 
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where  A ax  is  the  same  averaged  perturbation  equation  for  the  non-tilted  dipole  given  in  Eq.(75)  and 


2C 

A an  = asin;  j  7(^v)(l  +  2ecosv) 


2  nrj 


V  V  | 

cQcos — l-i’Qsin—  \dv 


i)  i j 

The  sine  and  cosine  functions  of  the  longitude  of  ascending  nodes  do  not  change  much  over  the  2 7tr] 
interval  and  may  be  considered  constants  for  the  integration  (they  are  abbreviated  as  sinQ  =  sQ.  and 
cosQ  =  cQ  for  clarity).  When  ij  is  an  integer,  the  only  terms  in  the  current  (Eq.(25))  that  will  resonate 
with  the  lower  frequency  harmonic  and  survive  the  integration  are  the  u6  and  u7  terms.  Therefore,  in  this 
case  the  contribution  to  the  change  in  the  average  semi-major  axis  due  to  the  rotating  tilted  dipole  is 

.2/Z77 


(26) 


A  an  =  CIma  sin  i  (cQm6  +  sClu1  j : 


(27) 


where  we  have  assumed  that  frequencies  of  the  F,  and  contributions  are  commensurate  and  that  terms 

with  incommensurate  frequencies  drop  out  after  integration  in  Eq.  (26)  yielding  an  exact  solution  shown  in 
Eq.  (27).  For  a  satellite  orbiting  at  an  altitude  of  approximately  261  km,  the  scaling  parameter  corresponds 
to  16  orbital  revolutions  per  sidereal  day,  i.e.  rj  =  16  .  To  consider  an  orbit  transfer  at  altitudes  that  do  not 
correspond  to  an  integer  number  of  revolutions  per  sidereal  day  where  the  multiple  frequencies  under 
consideration  are  not  commensurate,  we  average  the  state  over  a  larger  integer  number  of  revolutions  to 
achieve  an  approximate  model  for  the  averaged  state  dynamics.  This  is  accomplished  by  recognizing  that 
for  some  tolerance  r  >  0  ,  3  N  e  Z  and  p  e  Z  such  that  \i]N  -  p\<r  .  Simply  said,  if  we  choose  N  such 
that  an  interval  27rr]N  is  very  close  to  an  integer  number  of  periods,  then  the  commensurate  frequency 
model  in  Eq.  (27)  will  suffice  to  represent  the  averaged  dynamics  within  a  tolerance  that  is  defined  by  r  . 
This  means  that  the  duration  of  the  maneuver  must  be  long  enough  to  obtain  an  accurate  average  of 
instantaneous  states  that  include  contributions  at  lower  frequencies.  By  choosing  intervals  that  do  not 

correspond  to  integer  periods  (i.e.  rjN  <£  Z ),  the  maximum  mean  square  error  of  our  estimate  for 

a 

Cl 

incurred  by  using  Eq.(27)  is  of  the  same  order  as  the  nondimensional  quantity  — —  (~10'5  for  the  examples 

n 


here).  The  error  due  to  the  approximation  is  itself  periodic,  free  of  secular  growth,  and  is  exactly  zero 
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when  r/N  e  Z  .  Therefore,  after  integrating  Eq.  (26)  over  the  interval  from  0  to  2ktjN  ,  Eq.  (27)  becomes 
(for  N  »  1 ) 


=- 


ClaN 


2x1]  sin  i ( cQu6  +  sQu, )  + 


Similarly  the  average  inclination  change  due  to  the  contributions  of  the  rotating  tilted  dipole  is  derived 
using  the  method  of  averaging. 


2ni]N 


\  7(v') 


COS  cc 


cos  a  cos i 


'  V  V ' 

cQ  cos  —  +  sQ  sin  — 


-sin  a 


CT, 


2nr\N 


cos  i  cos  a 


v  .  v 
cd  cos  —  +  sQ  sin  — 


v 


V  .  V 

sQcos - cQsin  — 


\dv 


\dv 


Integrating  yields 


A-  CIm  (  m  m  \2xriN 

A i  « - cos  i  cilu,  +  silu1  - 

"  4  V  6  1J  n 


For  the  longitude  of  the  ascending  node,  the  change  due  to  this  contribution  is 


AQ„  = 


-  j  I{v)* 


2nrjN 


0 

2tujN 


A  j  si 

n  J 


sin  a 


r 


I  V  V 

cos  a  cos/ 1  cQcos  — +  sQsin  — 

> i  n) 


-sin  a 


f  v 

sQcos - cQsin  — 


dv 


2  I  V  .  V 

sin"  a  sQcos - cQsin  — 

V  ’1 

2m]N 


■  C 


V 

i  — 

7 


.  i^ 
L  — 

>U 


dv 


i  (sQu6  -  cQ.il-, ) 

4  n 


The  eccentricity  vector  components  are  derived  as  follows. 

1  'l  .  /  „  x  rcot /sin a  , 

slnr(Fli  et  j - 5 


2jrrjN  1 

f  J 

f* 

r 1  r 

J 

0  1 

l 

nae\ 

A  1  + 


ecosv 


- 

+h-C\ 

(  2  \  ] 

“ - r  |(F1  ’®t)  I 

na  e 1 

I  (l-2ecosv) 


2T-\,  ,  f  V 

k\ (2  +  ecosv)sinvsin/  cQ cos  —  +  sQ sin  — 

o  l  V  n  n) 


-sin  a 


I  V  Vi 

cos«cos/|  ciTcos  —  +  sQsin—  |  —  sin  a 
rj  77 


v  v 

sQ  cos  —  -  cQ  sin  — 


>I(v) 


+ /zsin/^2cosv  +  e(cos2  v  +  l)j  cQcos  —  +  s£2sin  — 
Integrating  yields 


Cl  ~ 

A hn  ~  ™  {Jc  cot  i ( sQu6  -  cQu 7 )  +  3 h  sin  i  (cQu6  +  sQu 


,  \  2  kijN 


n 


dv 
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2mjN 


Akn  = —  J  h j  (2  +  ecosv)sinvsin/ 


f  V  V  ^ 

cQ  cos  —  +  sCl  sin  — 


+coti  sin  a 


cosacos/j  cQcos  —  +  s£Tsin  —  |-sina 


^  V  V  ^ 

sCl  cos - cQsin  — 


-£sin/(2cosv  +  e(cos2  v  +  l)) 


>/(v) 


V  V 

cCl  cos  —  +  sCl  sin  — 

n  nj 


I (v)d\ 


which  produces 

Ak n  k,  ^^-^cot  i{-s£lu6  +cQi/7)  +  3A:sin/(cQM6  +  ,?Qm7  ) j 

For  a  given  averaged  state  x ,  the  total  rate  of  change  is  approximated  as 

Ax  «  Ax,  cos  S  +  Axn  sin  8 

where  the  non-tilted  dipole  dynamics  periodic  over  a  single  orbit,  Ax, ,  are  given  in  Appendix  G. 

2nnN  .  Ax 

Recognizing  that  AT  = - - —  ,  the  averaged  dynamics  -  may  be  determined.  For  a  sufficiently  long 

n  AT 

orbit  transfer  using  an  electrodynamic  tether,  this  averaging  method  will  capture  the  averaged  effects  due  to 
the  lower  frequency  rotation  of  the  tilted  dipole.  The  following  example  will  demonstrate  this  idea. 


Solution  to  an  Optimal  Control  Problem  Using  Multiple  Time 
Scales 

Using  the  tether  model  from  the  previous  optimal  maneuver  example,  a  longer  term  optimal  orbit 
change  maneuver  is  investigated  that  includes  the  moderately  varying  effect  of  the  Earth’s  rotating  tilted 
dipole.  For  this  model,  we  use  the  magnetic  field  described  by  Eq.  (23)  and  a  dipole  tilt  of  1 1.5  degrees. 
The  example  maneuver  will  increase  the  inclination  of  an  EDT  in  a  261  km  parking  orbit  (where  tj  «  16  ) 

from  40  to  45  degrees  ending  at  the  same  altitude  while  maintaining  a  constant  eccentricity  of  0.005  in  a 
drag  environment.  The  optimal  control  problem  is  therefore  written  as 
Minimize  Cost:  J  =  tf 

Subject  to: 
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x(r)  =  f(x(r),u(r)) 
e0  (x(r0 ))  =  [6639  km, 0,0.005, 40”, 0°f 
ef  (x(7}))  =  [af,ifJ  =  [6639  km,45°f 
g1(u(7’))  =  /l  -  2.25  <0  Amps2 
g2(x(r))  =  h2+k2-e20=  0 


where  the  mean  square  current  is  defined  using  Eq.  (13)  as 


rms  m 


2  1  /  2  2  2  2  2 

+  —^2  ^3  “I"  ^4  +  ^5  +  +  li 


v) 


The  dynamic  constraints  are  given  by  f  (x(r),u(r))  =  dx/dT  «Ax/A  T  and  box  constraints  given  in  Eq. 

(17)  are  also  enforced.  Solving  the  optimal  control  problem  using  DIDO  yields  the  optimal  control  profile 
shown  in  Figure  16. 


Min  Time  Orbit  Change  Control  Fourier  Coefficients  with  Tilted  Dipole 


Control  in  Time  Domain 
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Figure  16.  Control  Profile  Using  Tilted  Dipole  Model 

The  controller  dc  component,  w  ,  and  the  Fourier  coefficients  corresponding  with  the  higher  frequencies 
(  m,  through  u5 )  look  similar  to  the  corresponding  minimum  time  problem  in  the  previous  section  (see 
Figure  12).  Slightly  more  power  is  dedicated  in  the  form  of  direct  current,  corresponding  to  ul ,  to  change 
the  altitude  because  at  this  higher  inclination  orbit  the  local  magnetic  field  is  less  effective  for  thrust.  In  the 
tilted  dipole  case,  however,  a  small  controller  contribution  at  the  lower  frequency,  u6 ,  is  evident  which 
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superimposes  a  lower  frequency  component  on  to  the  control  signal  (see  bottom  of  Figure  16).  The  altitude 
and  inclination  trajectories  shown  in  Figure  17  reveal  a  similar  “climb  and  descend”  strategy  to  that  of  the 
example  in  the  previous  section.  This  maneuver  is  more  aggressive  than  the  previous  one  taking  113  days 
to  complete  and  close  inspection  of  the  propagated  trajectory  reveals  the  impact  of  a  rotating  tilted  Earth 
magnetic  dipole.  The  magnified  inserts  in  Figure  17  show  the  effects  of  the  3  time  scales;  the  fast  time 
dynamics  of  the  instantaneous  altitude,  the  medium  time  attitude  dynamics  with  daily  oscillations,  and  the 
slow  trend  of  the  average  altitude. 


State  trajectories  using  42  nodes 


Revolutions 

Figure  17.  State  Trajectories  Using  Tilted  Dipole  Model  with  Drag 

For  comparison,  optimal  controls  were  determined  for  the  same  problem  using  a  model  that  does 
not  include  a  tilted  magnetic  dipole.  The  controls  (Figure  18)  look  similar  to  the  previous  ones,  albeit 
without  the  medium  time  scale  components.  The  maneuver  appears  to  take  two  fewer  days  to  complete 
when  power  is  not  directed  to  compensate  for  the  magnetic  dipole  motion;  however  the  propagation  of  the 
altitude  does  not  match  very  closely  with  the  output  from  the  model  (Figure  19)  for  such  a  long  term 
maneuver.  The  propagation  was  performed  in  the  same  manner  as  the  previous  example  for  comparison, 
indicating  that  the  errors  are  model  errors  and  not  numerical  errors. 
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Summary 


Because  an  EDT  draws  its  thrusting  capability  from  the  Earth’s  magnetic  field,  it  is  important  to 
use  a  magnetic  field  model  with  an  appropriate  amount  of  fidelity.  Engineers  may  obtain  control  strategies 
using  a  less  accurate  model  for  maneuvers  that  do  not  span  many  revolutions  (in  this  example,  less  than 
about  300  revolutions),  however  for  transfers  that  take  a  very  long  time  a  tilted  dipole  model  must  be 


Min  Time  Orbit  Change  Control  Fourier  Coefficients  with  No  Tilted  Dipole 


Control  in  Time  Domain 


Figure  18.  Optimal  Control  with  Drag,  but  No  Earth  Dipole  Tilt 

considered.  Using  the  multiple  time  scale  approach  and  the  method  of  averaging,  one  can  include  this  low 
frequency  effect  in  the  model  by  introducing  a  new  time  scale  variable  in  the  controller.  The  next  step  to 
improving  the  controller  is  to  include  the  librational  motion  of  the  EDT  in  the  dynamic  model,  which  is  the 
subject  of  the  next  chapter. 
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Altitude  Trajectory  with  No  Medium  Time  Scale  Control 


Revolutions 


Figure  19.  The  controller  model  breaks  down  when  Earth  magnetic  dipole  tilt  is  excluded 
for  a  long  term  orbit  transfer.  Stars  indicate  the  model-derived  altitude  trajectory;  line  indicates 
propagated  altitude  trajectory  in  a  rotating  tilted  dipole  magnetic  field. 
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V.  Tether  Libration 


When  controlling  an  electrodynamic  tether  (EDT)  to  reach  a  new  orbit  as  discussed  in  Chapters  III  and 
IV,  it  was  assumed  that  the  tether  was  nadir-pointing  and  non-librating.  This  was  done  to  introduce  the 
method  of  averaging  for  solving  the  optimal  control  problem  in  Fourier  space.  In  reality,  however,  we 
would  need  to  account  for  the  librations  of  the  long  tether.  It  is  well  known  that  an  unperturbed  inert 
(unpowered)  tether  librates  in  and  out  of  plane  about  an  equilibrium  point  for  circular  orbits  without  growth 
or  decay.52'53  An  uncontrolled  EDT  with  a  constant  current  running  through  it,  however,  will  eventually  go 
unstable  as  aptly  pointed  out  in  Ref  10.  The  purpose  of  this  chapter  is  to  analyze  the  stability  of  an  EDT 
and  to  use  the  resulting  stability  criteria  to  define  libration  constraints.  The  objective  will  be  to  determine 
the  optimal  control  that  will  maneuver  an  EDT  to  a  new  orbit  while  simultaneously  driving  libration 
amplitude  to  a  desired  end  state  within  these  specified  libration  bounds.  To  achieve  this  objective,  this 
chapter  will  first  provide  an  examination  of  the  stability  of  the  tether  libration  both  with  and  without  drag. 
Given  that  an  EDT  with  a  dc  control  current  eventually  goes  unstable,  it  is  shown  that  the  system  may  be 
stabilized  using  a  method  of  feedback  linearization.  This  demonstration  provides  confidence  that  there  is  at 
least  one  feasible  control  solution,  thereby  allowing  us  to  seek  an  optimal  control  solution.  The  remainder 
of  the  chapter  is  devoted  to  the  derivation  of  dynamic  model  and  path  constraints  and  then  determining 
optimal  controls  for  a  librating  EDT  in  orbit  transfer. 

Equilibrium  and  Stability 

The  first  step  in  stability  analysis  is  to  obtain  the  attitude  dynamic  equations  of  motion  for  the  system. 
The  attitude  equations  of  motion  will  initially  be  based  on  the  following  assumptions  for  an  EDT  system  of 
two  subsatellites  connected  by  a  wire  in  tension.  These  assumptions  and  approximations  may  be  relaxed  as 
need  arises,  but  the  ones  listed  here  are  necessary  to  model  the  EDT  system  and  clearly  demonstrate  the 
utility  of  multiple  time  scale  optimal  controls  applied  to  libration  control. 
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Rigid  Rod  in  Tension  -  The  tether  is  assumed  to  be  perfectly  straight  between  two  subsatellites.  This 
approximation  is  valid  for  certain  ranges  of  maximum  allowable  wire  control  current  (see  Appendix 
D).  Because  the  tether  is  straight,  the  center  of  mass  (COM)  is  located  along  tether.  The  tether  cannot 
undergo  compression  or  go  slack,  but  rather  it  remains  in  tension  and  does  not  stretch.  The  former 
condition  is  valid  because  tether  attitudes  will  be  constrained  through  active  control  to  remain  below 
libration  angles  that  would  permit  a  slack  tether.  The  no-stretch  condition  is  justified  since  the 
materials  used  will  be  such  that  the  stretch  dynamics  is  insignificant  and  may  be  ignored. 

Medium  Length  Tether  -  The  tether  is  long  enough  to  consider  gravity  gradient  and  aerodynamic  torques 
due  to  air  density  variations  along  the  length  of  the  tether  to  be  significant.  The  latter  assumption  may 

be  restricted  to  W  *  «:  1  when  appropriate,  but  the  term  will  be  retained  for  generality  in  the 
/  h 

derivation  of  the  equations  of  motion.  The  characteristic  (or  scale)  height  of  the  atmosphere,  h  ,  is 
about  30  to  60  km  for  altitudes  between  150  and  400  km  [MSIS  Standard  Atmosphere].  See  Figure  20 
for  MSIS  standard  atmosphere  plots.  The  tether  is  considered  short  enough,  however,  the  magnetic 
field  is  approximately  constant  along  the  tether  length.  Implicit  in  this  assumption  is  the  tether  length 
is  small  compared  to  the  distance  to  the  center  of  the  Earth,  r ,  such  that  r  »  L  . 
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Maximum  Atmospheric  Density  Scale  Height 


Figure  20.  Standard  MSIS  Atmosphere 

Spherical  Earth  with  Non-tilted  Magnetic  Dipole  -  Although  the  magnetic  dipole  is  actually  tilted 
approximately  11.5  degrees  from  true  north  and  rotates  once  per  day,  this  effect  is  ignored  without 
severe  impact  to  the  initial  stability  analysis  and  control  design.  Figure  21  depicts  the  coordinates  used 
to  describe  the  in-plane  and  out-of-plane  librations  respectively. 
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Tether  Tip  Direction 


Orbit  path 


* 

Nadir  Direction 


Figure  21.  EDT  Attitude  Geometry  defining  the  in-plane  and  out-of-plane  libration  angles 

6  and  (j) ,  respectively 


It  may  be  desirable  to  have  the  tether  maintain  certain  attitudes  ( 9,  <!> )  or  operate  within  limits  of 

acceptable  attitudes.  With  the  equations  of  motion  we  can  proceed  to  determine  the  equilibrium  points, 
their  stability  and  the  non-linear  motion  of  the  tether  in  the  rotating  frame  of  reference.  The  libration 
equations  of  motion  were  derived  in  Appendix  A  using  Lagrange’s  method,  shown  here  employing  the 
rigid  tether  assumption. 


0  =  -t>'  +  2((9  +  v>)(zitan(zS-3:^f  sin  <9  cos  #  + 


Qe 


/ueL  cos  (/) 


(28) 


(j  =  -{(6  +  v)2  +  3^cos2  #}sin^cos(Z>  +  -^y 


(29) 


Variable  v  is  the  true  anomaly,  L  is  the  tether  length  and  ft,  is  the  effective  reduced  mass  (defined  in 

Appendix  B)  that  accounts  for  the  end-masses  and  the  tether  mass.  The  scalars  Qg  and  06  are  the 

generalized  forces  due  to  the  combination  of  electromagnetic  Lorenz  and  aerodynamic  drag  forces.  These 
equations  make  no  assumptions  about  the  ellipticity  of  the  orbit  and  may  be  related  to  the  rate  of  change  of 
the  true  anomaly  by 
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vg  _  V 

/'3  1  +  ecosv 

where  e  is  the  eccentricity  of  the  orbit.  We  start  by  analyzing  the  unperturbed  system  and  then  later  add 
some  of  the  perturbing  effects  like  atmospheric  drag. 


No  Drag  Model 

Unperturbed  tether  system  stability  has  been  analyzed  by  others9'10’52,  but  will  be  repeated  here  in 
a  manner  that  serves  the  purposes  of  this  research.  Starting  with  the  equations  of  motion  we  can  readily 

observe  the  equilibrium  points,  9e  and  <j>e ,  where  9  =  9  =  (/>  =  (/)  =  L  =  0. 


c  <f 


3  ju 

v  +  — f-c9s9 

\  r  j 


=  0 


C(j)S(l) 


(  3  u  ^ 

v2+-^fc29 

\  r  j 


=  0 


(30) 


(31) 


K 


Although  the  above  equations  allow  for  an  equilibrium  point  at^(,  =  ±—  when  the  tether  is  perpendicular  to 

the  orbital  plane,  we  shall  soon  discover  that  the  tether  cannot  maintain  pure  positive  tension  at  this  attitude 
thereby  allowing  the  end  bodies  to  orbit  separately  without  constraint.  In  controlling  space  tethers,  we  will 
avoid  this  case  since  we  desire  to  maintain  tether  tension  to  keep  a  valid  dynamic  model.  Other  equilibrium 
points  are  present  when  we  consider  a  circular  orbit.  With  a  circular  orbit,  the  true  anomaly  changes  at  a 

IjU 

constant  rate  i>  =  a>0  =  ,  —f  so  v  =  0  .  With  this  assumption,  consider  the  in-plane  libration  case,  <j>e  =  0  such 
V  r 

that  Eq.  (30)  reduces  to  3co^c9s9  =  0  while  the  second  is  satisfied  for  all  9 .  The  equilibrium  points  in  this 

7T  . 

case  occur  when  the  tether  is  in  a  lead-trail  co-orbital  configuration,  i.e.  9e  =  ±—  ,  or  in  a  nadir/zenith¬ 
pointing  configuration,  i.e.  (9e,<f>e )  =  (0,0)  or  (n, 0) .  It  will  be  shown  later  that  a  tether  can  go  slack  in  a 

lead-trail  orientation,  so  we  will  instead  avoid  this  configuration  and  only  investigate  the  system  stability  of 
the  nadir-pointing  equilibrium  point  where  positive  tension  can  be  maintained. 
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The  unperturbed  equations  of  motion  for  a  tether  in  a  circular  orbit  are  first  rearranged  as  explicit 
solutions  for  the  libration  accelerations. 

0  =  2  +  a>0  )  -  3  co^cdsO 

<j>  =  -c(j>s(j>  ^(<9  +  coa  )  +  3a>gC~0  j 

Defining  the  state  vector  as  x  =  [#,  <j>,  0,  we  can  generate  the  state  vector  time  derivative  and  its 
Jacobian. 

0 
</> 

l&ttj)  (  0  +  coo )  -  3  ofcBsB 
-c</)s</){^(0  +  a>oy  +  3o/c20j 

0  0  10 

0  0  0  1 

3  af(c20  -  sz0)  2  sec2  <)<)(0  +  a>o)  2^t(j)  2t<j>(8  +  wi) 

6 co^cBsBc^scj)  ~{c2<f>  -  s2(f>){(0  +  a>o)2  +3g)2oc10}  -2c(/>s(f>(8  +  a>o)  0 

Linearizing  A  and  evaluating  it  at  equilibrium  point  ( 0  ,  (/>  )  =  (0, 0)  yields 

0  0  10 

0  0  0  1 

-3®2  0  0  0 

0  -4®2  0  0 

Defining  a  =  -3 a>~  and  b  =  -4 co2  the  characteristic  equation  is 
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l4-(b  +  a)/l2+ab  =  0 


(34) 


|  U-A 


10  10 
010  1 
a  0 1  0 
Oft  01 


I2  =  b,  l2  =  a 

The  eigenvalues  are  therefore  ±yfb  and  ±Va  .  Note  that  as  long  as  a  <  0  and  b  <  0  then  the  system  will 
have  marginal  stability.  When  a  >  0  or  b  >  0  ,  then  there  will  be  a  positive  real  eigenvalue  which  indicates 

an  instability.  In  this  case,  the  eigenvalues  are  ±\3a>  i  and  ±2 co  i ,  pure  imaginary  numbers.  This  means 
that  in  the  vicinity  of  the  nadir-pointing  equilibrium  position,  the  tether  will  have  pendular  motion  with 
frequency  \[ia>o  in  the  orbital  plane  and  2a>  out  of  plane.  Only  circular  orbits  have  been  considered  at 

this  point.  For  non-circular  orbits  the  true  anomaly  rate  changes  with  respect  to  time,  therefore  the  system 
would  be  non-autonomous.  Floquet  theory  would  be  better  suited  to  determine  the  stability  of  this  system 
with  a  periodic  solution.  Palaez  et.  al  offer  a  more  thorough  discussion  of  the  stability  in  Refs  9  and  10  for 
a  powered  EDT  not  subject  to  drag. 


Drag  Model 

For  purposes  of  controlling  tether  libration  in  a  circular  orbit,  the  strategy  requiring  the  least 
amount  of  energy  would  be  one  that  controls  about  the  equilibrium  point.  When  atmospheric  drag  is 
considered,  the  equilibrium  point  may  be  slightly  different  than  that  of  the  tether  in  a  pure  vacuum.  To 
determine  this  equilibrium  point,  we  again  write  the  equations  of  motion  that  include  torque  due  to  the 
drag. 

nfic  <j>{0  +  3 w]cOsO  -  {0  +  <»  ))  =  Q0a 

juL2  +  cif>si/>((9  +  coo )'  +  3®2c26>))  =  Qia 

The  in-  and  out-of-plane  torques  for  a  tether  are  derived  in  Appendix  A  with  the  results  shown  here 
assuming  a  circular  orbit. 
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^  (p,  -1))} 

(35) 

(36) 

Qga  =  p{h)v2cc>ce^(By' -B'2ep>)-c(ep>  {-p2  -l)-e*  (Pl  -1)) 
Q»  =p(h)v2s<l>s0^(B;e» -B'2e-»)-C(e-»  (~Pl  -l)-e»  (p,  -1))} 


where  C  =  d  (l  -  cos2  ^sirf  (#  +  r)  p  2  -  and  the  non-dimensional  parameters  p,  and  p?  are 

given  as  p,  =  and  p ,  =  .  The  other  parameters  in  Eqs.  (35)  and  (36)  are  the  system 

/?  Mj  /?  I/  , 

velocity  V,  the  atmospheric  density  p  (  h  )  at  altitude  h  ,  and  B  values  representing  the  ballistic 


coefficients  of  the  end-masses.  The  mass  parameters  M  and  M,  are  defined  as  M  =  in  h — and 

i  i  2 


M  =  r,  h — 1  where  ni  is  the  mass  of  end-body  1,  »i,  is  the  mass  of  end-body  2,  and  w  is  the  mass  of  the 
2 

tether  (see  Appendix  B  for  details). 

To  obtain  equilibria  of  long  tethers  at  various  altitudes,  no  simplifying  assumption  on  the  size  of 

W ,  has  been  made  yet.  According  to  the  equations  of  motion  equilibrium  is  achieved  when  (/)  =  (>  =  0 

/  h 

which  occurs  when  <j>  =  0  and  6  satisfies  the  following  equation. 

pmLcde  -pmLc9e  \  Z'  - pm Lc9e  /  x  jumLc9e  /  \  \  'j 

P„Z  n*  M  h  n*  M  h  M  h  ~ Pm  LC’9  M  h-  [A  Lc9  | 

-  B  e  -  Bne  -Ce - 1  -e - 1  f 

2  '  Mh’  Mh’ 


M,  h  ym,  M,h  Pm  Lc9  I  < 

-  C  e - 1  —e 


p  Lc9 


Recall  that  B  =  1/1  1  and  B ,  =  Ji  2  and  when  </>  =  0  in  a  circular  orbit,  then  C  =  d  —  .  This  result  is 

Ml  M2  ‘  c9 

consistent  with  that  of  Beletsky  and  Levin  [Ref  52,  p  214  and  262]. 

This  indicates  that  an  equilibrium  point  resides  in  the  plane  of  the  circular  orbit  and  is  offset  from 
nadir  pointing  by  an  angle  9  that  satisfies  the  transcendental  Equation  (37).  Solutions  to  Eq.  (37)  are 
obtained  numerically  for  given  values  of  altitude,  density  and  tether  characteristics.  Figure  22  and  Figure 
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23  show  the  equilibrium  points  residing  in  the  orbital  plane  for  various  tether  characteristics  and  altitudes. 

These  plots  show  three  main  trends.  They  are 

•  Increasing  the  altitude  drives  the  equilibrium  point  to  nadir.  In  the  absence  of  drag,  the 
equilibrium  point  is  exactly  zero. 

•  Increasing  the  disparity  between  the  upper  and  lower  endmasses  drives  the  equilibrium 
point  away  from  nadir.  When  one  mass  is  more  massive  than  the  other,  it  is  less 
susceptible  to  drag  resulting  in  a  more  tipped  orientation  on  average. 

•  Increasing  the  tether  length  beyond  3  or  4  km  up  to  10  km  does  not  significantly  affect 
the  equilibrium  point. 

L 

It  we  make  the  approximation  that  «  1 ,  then  the  equilibrium  condition  reduces  to 

h 


In  this  form,  we  may  explicitly  solve  for  9e  .  Caution  must  be  exercised  since  s9e  is  bounded  by  [-1.1] 

(*  *  \  . 

B{  -B2  1  is  also  bounded  to  certain  values  for  given  tether 

dimensions  and  mass  distribution.  There  is  a  nadir-pointing  equilibrium  condition  where  (Qe,$e )  =  (0, 0) 

* 

(*  *  \  dth 

5,  -  B2  I  »  -2  — —  . 

Mm 

Stability  about  the  neighborhood  of  an  arbitrary  equilibrium  point  may  be  determined 
approximately  using  the  same  technique  as  was  done  in  the  case  with  no  drag.  When  the  in-plane  angular 
acceleration  due  to  aerodynamic  drag  torque  (Eq.  (62)  in  Appendix  A)  is  included  in  the  dynamics,  the  state 
vector’s  time  derivative  is  given  by 
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f  =  i  = 


2<j>t$[6  +  <voJ-3tv2c9s9  +  — ^  ^ |^- a/)c9 (fi*  -  B^  +  dJi  (\-c2(/)s20y 
-c(/)s(/>  ((<9  +  co0  )2  +  3®;c26>  j  +  V  scj)s9  \^~  (#*  -fi,*)- (l  -  c  V-r  6>) 


In  calculating  the  linearized  A  matrix  about  the  in-plane  equilibrium  point  (4),  0)  with  0e  =  (f>e  =  0 ,  the 
only  terms  that  will  differ  from  those  given  in  Eq.  (33)  are  An  =  and  A4 ,  =  .  Because  of 


the  structure  of  the  linearized  A  matrix,  the  eigenvalues  will  have  positive  and  negative  real  parts  if  these 
two  elements  of  the  matrix  are  positive  and  an  instability  would  exist  (see  the  characteristic  Eq.  (34)).  If 
these  elements  remain  negative,  then  there  would  be  a  pendular  libration  in  the  neighborhood  of  the 
equilibrium  point,  as  we  had  with  the  no-drag  case.  The  first  element  under  investigation  is  given  as 


4, 


We  may  substitute  the  equilibrium  condition  (38)  into  the  second  term  on  the  right  hand  side  of  A}1  to 


produce 


4,  =  -34 (c29e  - s%)-  =  -34 (c29e  - s2de) - 3 co2s29e 


An  =  -?>(D-^6e 


Likewise,  the  A  term  is  calculated  as 


5/4  =  -Acolc2d{c l<f  - 5 V)  +  v~P(h)c*s0  W  (B*  _  B- )  _  (1  _  c  V*)' 

d<j>  <ueL  \  2  x  '  cdc^  ’ 


,  2  2.  v  Pi 
=  -4®n  C  ~ - 


p(4  l  2 


Once  again,  substituting  the  equilibrium  condition  (38)  into  the  second  term  on  the  right  side  yields 


4,2  =  -4colt26e  - 3a>*s20e 


4,2  -  w0 


(3  +  c2ee) 
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Notice  that  both  terms  A3I  and  A42  remain  negative,  thus  in  the  vicinity  of  the  equilibrium  point  the  tether 

librates  with  in-plane  frequency  ^3 cOecoo  and  out-of-plane  frequency  ^3  +  c~9ea>0 .  For  the  nadir-pointing 
case,  we  have  the  same  solution  as  that  of  the  no-drag  case. 
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Tether  Length  [km]  Tether  Length  [km] 


In-plane  Lib  angle  equilibrium  for  various  tether  designs,  180  km  alt  circ  orbit 


Equilibrium  Theta  [deg] 


In-plane  Lib  angle  equilibrium  for  various  tether  designs,  200  km  alt  circ  orbit 


Equilibrium  Theta  [deg] 


Figure  22.  In-plane  equilibrium  points  for  180  km  and  200  km  circular  orbits 
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Tether  Length  [km]  Tether  Length  [km] 


In-plane  Lib  angle  equilibrium  for  various  tether  designs,  250  km  alt  circ  orbit 


Equilibrium  Theta  [deg] 


In-plane  Lib  angle  equilibrium  for  various  tether  designs,  300  km  alt  circ  orbit 


Equilibrium  Theta  [deg] 


Figure  23.  In-plane  equilibrium  points  for  250  km  and  300  km  circular  orbits 
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Tension 


When  determining  the  equilibrium  point,  we  assumed  that  the  tether  was  rigid,  which  is  only  a 
good  approximation  when  the  tether  is  in  tension.  In  the  case  of  a  slack  tether,  we  would  have  to 
unconstrain  the  equations  of  motion  and  retain  the  L  term  in  the  equations  of  motion.  We  can  use  the 
equations  of  motion  to  determine  if  the  tether  is  in  tension  at  the  nadir  pointing  equilibrium  point.  The 
tension  is  depicted  in  Figure  24. 


Figure  24.  Tether  Tension  Diagram 


The  generalized  force  along  the  tether  length  is  derived  as  follows. 


Qu  =  T 
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fri 

dL 


5v,  5v,  -ju„,  dv, 

-  t - —  where  t  =  rp,  — —  = - p,  and  — — 
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T  =  ~T 
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The  force  due  to  the  atmospheric  drag  along  the  tether  for  (j)  =  0  is 


{h}v2sO^ 


BxePx  -B2e  Pl 
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The  electrodynamic  force  does  not  add  any  force  along  the  tether  since  it  acts  perpendicular  to  the  assumed 


Lu„ 

straight  tether.  So  from  the  equation  of  motion  in  L  for  a  circular  orbit  (Eq.  (59)),  i.e.  v  =  =  a>o ,  and 

r 


rigid  rod  (Z  =  Z  =  0)  in  equilibrium  we  have 


ml 


L 


f  - c-<t> ( 0  +  v)2  -  3c2<tc20  - 1)  =  -3 =  QLa  +  Qu 

r 


(39) 


-3jueLa?c28  =  ^'L p(h)v2s0^BlePl  - B2e  P2^j-T 

r  =  3jueLo)2oc20  +  ^-p{h)v2sd(yBlePl - B*2e~Pl )  (40) 

This  result  is  consistent  with  Ref  45  p.  125.  We  defined  r  to  be  in  tension  when  it  is  positive  as  defined  in 
Figure  24.  The  first  term  on  the  right  hand  side  of  Eq.  (40)  is  always  positive  and  dominates  the  second 
term  near  9  =  0,  n  except  when  the  term  in  brackets  is  a  very  large  value.  This  would  occur  in  extreme 
cases  where  the  design  of  the  tether  is  such  that  the  upper  and  lower  masses  have  very  different  ballistic 
coefficients  so  that  one  end  mass  is  subject  to  much  more  drag  than  the  other.  For  example,  if  the  in-plane 
libration  angle  0  is  positive  (i.e.  the  upper  mass  leads  the  lower  one)  and  the  upper  mass  undergoes  so 
much  drag  that  it  falls  behind  faster  than  the  gravity  gradient  can  correct  to  the  vertical  position,  then  the 
tether  could  go  slack.  The  graphs  shown  in  Appendix  E  present  numerical  solutions  for  tension  T  at 

different  in-plane  equilibrium  points.  The  remaining  graphs  in  Appendix  E  are  9e  solutions  for  the  zero 

tension  condition.  Each  of  these  angles  will  serve  as  an  upper  bound  for  in-plane  libration  so  the  tether 
does  not  go  slack.  For  this  design,  the  tether  would  need  to  be  fairly  close  to  90  degrees  (i.e.  lead-trail 
configuration)  for  reasonable  ballistic  coefficients  before  losing  tension.  The  tension  in  the  nadir-pointing 

position  is  r  =  3p  Loo  .  For  a  2  km  tether  in  a  250  km  orbit,  this  force  would  be  about  0.7  N. 

The  other  singularity  points  mentioned  were  ruled  out  because  the  tether  cannot  maintain  positive 

K 

tension  in  those  circumstances.  For  the  singularity  at  </>-± —  we  have  the  equation  of  motion 
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r  =  "  — ^—  +  QLa  with  0La  =  0. 
r 

This  indicates  that  the  force  along  the  tether  wire  would  be  in  compression,  which  of  course  is  physically 

7T 

impossible.  For  the  circular  orbit  singularity  at  9  =±—  (lead-trail  orientation) 

2 

t  =  -jUeLa^s2tf>  ±  v2/?(/?)c(zi  (b*  -  B*2  ) 

Tension  in  this  case  is  only  positive  when  t/>e  =  0  (or  n),  6e  =  +  ^or  and  ^  >  52  ,  i.e.  the 

trailing  end-mass  in  the  lead-trail  formation  is  subject  to  a  greater  drag  force.  Otherwise  the  tether  goes 
slack.  The  stability  of  EDTs  has  been  explored  by  other  researchers  who  conclude  that  when  driving  a 
constant  uncontrolled  current  through  the  tether  wire,  the  system  will  eventually  go  unstable,  tumbling  end 
over  end.  Instead  of  reproducing  the  results  here,  we  will  explore  an  example  of  a  controller  that  uses 
feedback  linearization  to  drive  down  the  libration.  After  gaining  confidence  that  the  system  may  be 
stabilized  by  employing  active  control,  we  will  turn  to  posing  and  solving  optimal  control  of  a  librating 
EDT. 

Demonstration  of  Attitude  Control  Using  Feedback  Linearization 

With  the  dynamic  models  presented  in  Chapters  III  and  IV,  we  are  positioned  to  explore 
electrodynamic  tether  libration  control  strategies.  Before  determining  optimal  control  for  a  librating  EDT, 
in  this  section  we  will  first  determine  a  feasible  solution. 

Libration  control  example  using  feedback  linearization 

There  are  two  types  of  2-ball  tether  system  attitude  control  strategies,  a  hanging  tether  and  a 
spinning  tether.  A  hanging  system  will  use  active  control  to  maintain  pendular  motion  about  equilibrium 
whereas  a  spinning  system  will  allow  the  system  to  tumble  end-over-end,  thus  avoiding  the  need  for  active 
attitude  control.  In  this  section  we  explore  a  possible  attitude  control  strategy  using  torque  resulting  from 
the  applied  Lorenz  force.  Although  not  optimal,  feedback  linearization  provides  a  quick  way  to  see  if 
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attitude  control  is  achievable  for  a  given  system  using  nothing  but  the  Lorenz  force  on  the  wire.  We  start 
with  the  dynamic  equations  of  motion  derived  in  Appendix  A.  They  are  repeated  here  for  the  state  vector 


f  =  x  = 


e 

</> 


2  </>t<f>(d  +  co0 )  -  2>arc6sQ  - 


Qbi i  +  Qee 
,2„2a\  Qpa  Qpe 


+  3 a  c  0\  + 


A/ 


We  will  take  advantage  of  the  assumption  that  the  tether  length  is  significantly  shorter  than  the  atmospheric 

scale  height,  «  1  and  write  the  generalized  forces  due  to  drag  from  Eq.  (62)  and  Eq.  (63)  as  the 
h 

following. 


Qe„  =  v2yo(/?)L  ^-^-afcO^B*  -  B* )  +  dth*  (l-c2  </>s2  6>j'' 


iX  h 


Q*a  =v-p(h)Ls</>sO^{B‘2-Bl)-d,{\-c1  (fis1 6)  ^ 


The  generalized  forces  due  to  the  Lorenz  force  acting  perpendicular  to  the  tether  are  given  by  Eq.  (64)  and 
Eq.  (65)  and  repeated  here. 


Qric 

Q<he 


IL2  (m0  -  ml ) 
2  M 

IL2  (  m1  -  ml ) 
2  M 


( c(j>c6s(fBr  +  s(j>c(j)s6Bl 
(- sGBr+cGB, ) 


-c^Bn) 


Suppose  we  wish  to  minimize  the  error  2-norm  of  the  tether  attitude  with  the  equilibrium  attitude,  i.e.  we 
desire  to  drive  the  states  y  =  h(x)  =  to  some  position  close  to  the  point  \0e,</>e ]r  .  Looking  at  the 

dynamics  of  the  system,  the  control,  I  =  u(t) ,  appears  only  in  the  second  derivative  of  y  .  Therefore, 
following  the  discussion  in  Ref  49  pp.  267-232,  we  may  obtain  the  dynamic  inversion  by  taking  the  second 
Lie  derivative  of  h(x) . 
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Now  to  derive  the  controller,  we  write 


y  = 


y  = 


d\ 

dx 2 


=  ^Tf  +  G  (x)u 
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Let  the  desired  values  of  6  and  (f)  be  given  by  the  vector  v .  The  controller  will  drive  the  states  toward 


the  desired  values  using  a  gain  matrix  K  . 
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Where  the  error  vector  and  gain  matrix  are  respectively, 


yi~ydi 

y2-ydi 

Ti-Trfi 


and  K  = 


kx  0  k2  0 
0  k}  0  k4 


\_y 2  To'2  J 


Setting  the  desired  output  equations  containing  the  gain  equal  to  the  system  dynamic  equations,  we  can 
solve  for  u(t ) . 


52h 

y-  =v  =  T-rf+M 

ox 


u  =  g;\x) 


d~h  f 

v - -f 

dx2 


(41) 
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The  term  Gs  1  (x)  is  the  pseudo  inverse  of  input  vector/matrix  G  and  is  required  for  dynamic  inversion 
since  there  are  two  controlled  states,  but  only  one  input  control. 

Although  the  pseudo  inverse  cannot  drive  the  tether  attitude  to  the  exact  desired  equilibrium  point 
in  general,  it  does  minimize  the  error  2-norm.  The  proof  is  shown  as  follows. 

.  We  can  write  the  error  norm  as  e  =  ||Gm  -  q||‘  and  minimize  this  with  respect  to  the 

control. 


Let  q  = 


dhr 

v - -f 

5x 


f=±([G„-qnG„-q])=° 

=  [Gw  -  qf  G  +  Gr  [Gw  -  q]  =  0 
=  2GrGz/  -qrG  -Grq  =  0 

Since  q  and  G  are  both  vectors, 

2GrGz/-2Grq  =  0 

=>  «  =  (GrG)"1Grq  =  -^rq  =  G;1q 

INI 

The  control  law  given  by  Eq.  (41)  was  implemented  using  Simulink  and  the  system  was  given  a  small 
perturbation  from  equilibrium.  The  output  plots  in  Figure  25  show  that  the  controller  drives  the  in-  and 
out-of-plane  libration  angles  back  toward  the  equilibrium  point.  Although  the  controller  does  not  minimize 
energy  or  time,  it  does  demonstrate  that  there  is  potential  for  controlling  attitude,  if  required,  using  only  the 
current  in  an  electrodynamic  tether.  In  reality,  consideration  must  be  given  to  the  orbital  change  impacts  of 
attitude  controlling  in  this  way.  Other  researchers  have  investigated  tether  length  or  tension  control  to 
dampen  librations. 50  The  next  section  will  demonstrate  how  one  may  use  optimal  control  methods  over 
large  time  scales  to  maneuver  the  tether  system  to  a  new  orbit  while  constraining  the  libration  angles  and 
rates  to  desired  values. 


58 


Current,  I  [Amps] 


5 


Feedback  Linearization:  4  km  EOT  in  160  km  circular  10  deg  inclination  orbit. 
Tether  properties:  M=400  kg,  m1/M=.5,  m2/M=.3,  mt/M=.2,  dt=5  mm. 
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Figure  25.  EDT  Attitude  Control  Using  Feedback  Linearization 
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Libration  Control  over  a  Long  Time  Scale 


When  controlling  an  electrodynamic  tether  to  reach  a  new  orbit  as  discussed  in  the  Chapter  III,  it 
was  assumed  that  the  tether  was  nadir-pointing  and  non-librating.  This  was  done  to  introduce  the  method 
of  averaging  for  solving  the  optimal  control  problem  in  Fourier  space.  In  reality,  however,  we  would  need 
to  account  for  the  librations  of  the  long  tether.  It  has  been  shown  that  an  unperturbed  inert  (unpowered) 
tether  in  a  circular  orbit  librates  in  and  out  of  plane  about  an  equilibrium  point  without  growth  or  decay. 
An  uncontrolled  EDT  with  a  constant  direct  current  running  through  it,  however,  will  eventually  go 
unstable  as  aptly  pointed  out  by  Pelaez,  Lorenzini,  Lopez-Rebollai,  and  Ruiz  in  Ref  10.  The  purpose  of 
this  section  is  to  incorporate  constraints  on  libration  in  the  optimal  control  problem  of  Chapter  III  that  will 
enable  an  optimal  orbital  change  maneuver  while  simultaneously  driving  libration  amplitude  to  a  desired 
end  state  within  specified  bounds.  Unfortunately,  straightforward  averaging  of  the  derivative  of  the 
libration  angle  as  we  did  with  the  orbital  state  derivative  would  yield  zero.  Control  cannot  be  achieved  for 
a  state  that  is  always  zero,  so  a  different  approach  is  required  to  capture  the  librational  motion  in  Fourier 
space  to  control  the  averaged  state. 

To  simplify  the  problem,  in-plane  libration  is  ignored  and  attention  is  placed  on  controlling  the 
out-of-plane  libration.  In-plane  libration  is  not  resonant  with  the  periodic  controller  or  the  orbital  motion 
(recall  from  the  Equilibrium  and  Stability  Section  that  coe  =  \[3n  where  n  is  the  mean  motion  of  the 
satellite),  thus  it  does  not  grow  very  quickly.  For  the  design  proposed  here,  months  of  constant  thrusting 
are  required  to  gain  a  few  degrees  of  in-plane  libration  amplitude.  Furthermore,  the  small  in-plane 
librations  may  be  managed  by  other  mechanisms,  such  as  controlling  the  drag  on  the  upper  and  lower 
bodies  thus  imparting  an  aerodynamic  torque  out  of  phase  with  the  pendular  motion  thus  dampening  this 
motion.  With  this  justification  in  mind,  we  derive  a  new  state  that  captures  only  the  out-of-plane  libration 
(hereafter  simply  called  “libration”  unless  otherwise  stated). 

A  constraint  in  Fourier  space  must  not  contain  any  functions  of  a  fast  time  variable,  i.e. 
trigonometric  functions  of  v  .  Averaging  serves  to  eliminate  dependence  on  this  fast  time  variable  leaving 
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only  variables  changing  slowly  with  time.  To  accomplish  this,  a  new  state  is  devised;  the  mean  square 
value  of  a  tether’s  out-of-plane  libration.  Whether  power  is  applied  to  the  tether  or  not,  the  libration  mean 
square  is  proportional  to  the  maximum  angle  reached  throughout  the  pendular  cycle.  For  an  unpowered 
(inert)  librating  EDT,  the  mean  square  value  is  exactly  half  of  the  square  of  the  libration  magnitude,  i.e. 

2 

(jrms  =  .  This  relationship  is  approximate  for  a  powered  EDT  as  long  as  the  perturbation  due  to  the 

electromagnetic  torque  is  relatively  small.  Deriving  an  expression  that  describes  the  librational  mean 
square  behavior  provides  a  way  to  understand  the  behavior  of  the  magnitude  of  the  librational  motion  over 
a  long  time.  Thus  constraining  the  mean  square  trajectory  for  a  given  orbital  maneuver  is  tantamount  to 
bounding  the  envelope  that  contains  the  librational  motion  of  the  tether  over  long  time  durations. 

Unfortunately,  the  librational  equations  of  motion  given  in  Eqs.  (28)  and  (29)  have  no  closed-form 
solution  that  will  enable  us  to  capture  the  libration  amplitude  changes  over  long  time  scales.  The  good 
news,  however,  is  that  assuming  small  libration  angles  we  may  linearize  the  equations  of  motion,  thus 
decoupling  the  in-  and  out-of-plane  libration  equations  of  motion  and,  as  previously  mentioned,  ignore  the 
in-plane  libration.  Introducing  the  mean  anomaly  v  as  the  independent  time  variable  we  write 

<j>  +  4<j)  =  - i  )'/m  i m  sjn  i  fk(T)cosv  -  /7  (r)  sin  v  Wr  (v)u(r) 

MeL  2M  juefug  (42) 

=  s  sin  i  [k  cos  v-h  sin  v)  (w,  +  u2  cos  v  +  u2  sin  v  +  u4  cos  2v  +  u5  sin  2v) 


where  dots  indicate  differentiation  with  respect  to  v  ,  i.e.  (  )  = 


dv 


This  equation  is  expressed  using  a 


partial  equinoctial  element  set  described  in  Chapter  III  where  k  =  —  and  h  =  —  .  Both  k  and  h  are  order 

e  e 

one  quantities  that  are  themselves  averages  that  vary  slowly  with  time.  Adopting  the  convention  used  in 

Refs  9  and  10,  the  non-dimensional  small  parameter  e  is  defined  as  the  ratio  of  the  maximum 

electrodynamic  torque  to  the  gravity  gradient  torque  and  corresponds  to  the  powered  part  of  the  expansion. 

^  _  Max  Electrodynamic  Torque  _  ( m2  -  mx )  ym  ^ 

Gravity  Gradient  Torque  2 M jueju 
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For  a  1.5  Amp,  500  kg  tether  system  in  low  Earth  orbit  with  an  upper  end  mass  of  230  kg  and  a  lower  end 
mass  of  220  kg,  this  parameter  is  about  0.0026.  The  current  control  introduced  in  Chapter  III 

as/(u(T),v) ,  has  a  slowly  varying  part,  u(T)  =  [ul,u2,ui,u4,u5]'  and  a  periodic  part  that  forms  the  basis 

in  F ourier  space  'F (V)  =  fl, cos  v, sin  v, cos  2v. sin 2v]r  .  The  normalized  control  current  is  therefore  given 

by  /  =  ImW7  (v)u(J)  where  Tr(v)u(r)  is  an  order  one  quantity.  Recall  that  the  slow  time  scale  variable 

T  is  a  scaled  version  the  clock  time  t  and  the  true  anomaly.  It  is  now  necessary  to  formalize  the 
relationship  between  the  two  time  scales  using  a  scaling  parameter  e  such  that 

T  =  st  =  £—  (44) 

n 

The  non-dimensional  scale  factor  used  here  is  the  torque  ratio  defined  in  Eq.  (43)  (see  Appendix  F  for 
details  on  scaling).  Only  small  changes  to  the  known  periodic  libration  motion  of  the  inert  tether  over  short 
time  spans  will  occur  as  long  as  the  electrodynamic  torque  is  small  compared  to  the  gravity  gradient,  i.e. 

£  «  1 .  In  transforming  the  controls  from  the  short  time  scale  domain  to  Fourier  space,  we  exchange  a 
single  control  variable  (current  as  a  function  of  a  fast  time  variable)  for  five  control  variables  (the  five 
Fourier  coefficients  in  this  case  that  are  functions  of  slow  time  variables).  Expanding  the  right-hand  side 
term  in  the  differential  equation  in  Eq.  (42)  and  through  liberal  use  of  trigonometric  identities,  we 
determine  an  exact  solution  to  the  linearized  equation  applying  the  method  of  undetermined  coefficients. 

<j>(y,T)  =  ^>0(y,T)  +  £<j>x  [v,  u  (71))  =  </,„  cos  2  (v  -  )  +  £  sin  i  (kK  -  /?//)  (45) 


where 


A:Mr))  =  |+^  +  y) 


W,  .  U,  U-,  .  u.  u.  . 

cos  v  +  —  sin  v  -  —  v  cos  2v  +  —  v  sin  2v  -  —  cos  3v  -  —  sin  3 v 


10 


10 


H (l/jUfr))  =  —  +  —  cosv-—  —  -u,  sinv-  — vcos2v  +  — vsin2v  +  — cos3v-  — sin3v 
1  V  ”  8  6  31  2  1 J  8  8  10  10 


One  restriction  due  to  the  linearization  is  that  the  second  term  on  the  right  hand  side  of  Eq.  (45) 
must  be  less  than  order  one,  i.e.  £V  <3C  1 .  Therefore,  to  ensure  accuracy  of  the  solution,  the  duration  is 


limited  to  v  «  —  (note  the  explicit  v  terms  present  in  K  and  H  ).  For  a  scaled  maximum  electrodynamic 

£ 
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torque  of  e  =  0.0026  ,  this  maximum  allowable  duration  corresponds  to  about  60  orbital  revolutions. 
Eventually  a  long  duration  optimal  control  problem  will  be  discretized  into  smaller  intervals  that  are  much 
shorter  than  this  limit  so  that  the  approximate  solution  is  valid  for  each  subinterval.  Linking  the 
sub  intervals  together,  the  long  term  maneuver  consists  of  the  slowly  varying  “constant”  states  and  Fourier 
coefficient  controls  within  each  subinterval.  The  first  term  on  the  right  side  of  Eq.  (45)  represents  the 
homogeneous  solution  indicating  that  a  tether  without  any  electrodynamic  torque  would  continually  librate 
at  twice  the  orbital  frequency.  Perturbations  come  through  the  small  electrodynamic  torque  of  order  e 
imparted  on  the  tether  over  a  long  time.  Whether  these  perturbations  destabilize  or  stabilize  the  libration 
depends  on  the  slowly  changing  control  terms  contained  in  K  and  H  .  A  thorough  derivation  of  periodic 
libration  angle  solutions  is  provided  in  Ref  10  for  an  EDT  with  a  steady  dc  current. 

For  an  unpowered  tether,  or  one  where  the  center  of  mass  is  collocated  with  the  center  of  force  on 
the  tether  (thus  no  Lorenz  torque),  £U  (J)  =  0  ,  or  an  equatorial  orbit  where  i  =  0 ,  the  solution  to  Eq.  (45) 
is  the  homogeneous  solution. 

Qo  {V’T)  =  <!>„,  (7’)cos(2(v-v0)) 

where  (j)m  ( T )  is  the  initial  amplitude  of  the  librational  motion  which  is  approximately  constant  over  a 
period,  but  changes  slowly  over  time.  Presuming  that  the  periodic  control  may  be  started  at  any  time 
during  the  libration  cycle,  for  purposes  here  we  assume  the  peak  of  a  libration  cycle  corresponds 
with  v0  =  0  and  write 

0o(v’T)  =  </'m{r)cos2v  (46) 

Using  this  model,  the  only  way  to  control  the  libration  is  through  the  0(s )  term  in  Eq.  (45). 

We  can  now  define  the  libration  mean  squared  state  as 

i  v+2/r 

z(yJ)=2 -  J  *2(£>r)  </£  =  £,  (47) 

This  state  is  always  positive  and  is  itself  an  average  over  a  period  by  definition.  Furthermore,  for  short 
time  intervals  such  as  a  few  periods  when  the  libration  amplitude  change  is  negligible,  the  relationship 
between  the  state  z  and  the  amplitude  may  be  expressed  as 
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The  first  integral  term  in  Eq.  (50)  is  the  inert  tether  libration  mean  square  value.  Secular  changes  enter  the 
system  through  the  remaining  terms  with  explicit  dependence  on  v  .  Over  a  single  period  the  change  in 
z  is  very  small  due  to  the  scaling  factor  s  .  Recalling  Eq.  (44),  we  substitute  the  slow  time  variable  for 
ev  ,  consider  it  constant  over  the  limits  of  the  definite  integral,  and  remove  it  from  the  integrand.  This 
slow  time  variable  affects  the  secular  growth  (or  decay)  of  the  z  state  only  over  large  spans  of  time,  so  only 
the  sinusoidal  functions  of  v  are  averaged  through  integration.  Physically,  the  mean  squared  value  of 
libration  changes  approximately  linearly  with  T  to  first  order  over  short  time  intervals.  The  plot  in  Figure 
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26  depicts  the  small,  nearly  linear  change  in  the  z  state  over  one  period.  Substituting  the  slow  time 
variable  into  Eq.  (50)  and  expanding  yields 


1  ^  yiT  1  ^  / 

? (r)  =  —  j  <fldv  +  —  <j)m  sin i  —  j  ^/?z^2  -  ku3  j cos2  2v  +  [ku2 
J-U  0  4  2  n  0  v 


«2r2sin2/  1 
64  2k 


2v)  +  /?"  (iq  cos2  2 v  +  iq  sin2  2v)  -2 hku2u2  \dv  +  0[sT ' 


Finally,  we  perform  the  integration  with  respect  to  true  anomaly  and  take  the  derivative  with  respect  to 
clock  time  for  the  desired  secular  change  in  z  over  a  long  time  scale.  Assuming  that  the  averaged  states  x 
and  control  coefficients  u  in  the  previous  equation  vary  slowly,  the  derivatives  with  respect  to  clock  time 
will  be  small,  i.e.  dxjdt  =  0(xs)  and  dujdt  =  0{us)  .  The  z  state  derivative  is  therefore 


dz  dz  ~j2zn  sin  i 

dt  dT  8 


(hu2  -  ku: )  - 


2  •  2  • 
2  n  sin  i 
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((ul  +  Mj )  —  4 hku2u3  j t  +  O(0m£2t^ 


(51) 


where  T  =  st .  Although  the  second  term  on  the  right  hand  side  of  the  derivative  in  Eq.  (51)  causes 
quadratic  growth  (or  decay)  of  the  z  state,  it  is  of  order  s1  and  may  only  be  significant  when  considering 
larger  time  spans.  This  derivative  will  serve  as  a  dynamic  constraint  in  subsequent  optimal  control 
problems  to  manage  the  magnitude  of  libration  while  performing  orbital  maneuvers.  Notice  that  using  this 
model,  the  change  in  the  libration  mean  square  state  is  achieved  primarily  through  the  u2  and  u2 
coefficients  corresponding  to  periodic  control  resonant  with  the  orbital  frequency.  This  is  because  in  the 
satellite  frame  the  local  magnetic  field  vector  varies  with  the  orbital  frequency,  therefore  resonating  control 
current  with  this  frequency  can  dampen  (or  excite)  libration. 
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Libration  Squared  Function  and  Envelope 


Figure  26.  Libration  Squared  Function  and  Envelope 


Optimal  Maneuvering  with  Libration  Control 

With  the  dynamics  of  the  mean  square  libration  in  hand,  it  is  possible  to  optimally  maneuver  an 
EDT  satellite  to  a  new  desired  orbit  while  controlling  the  out-of-plane  libration  (within  the  limits  of  the 
dynamic  model).  The  assumptions  made  for  this  section  are  that  the  in-plane  libration  is  controlled  using  a 
separate  mechanism  (e.g.  drag  torque  control)  and  that  the  out-of-plane  librations  are  much  larger  than  the 
in-plane  librations,  i.e.  0  <^<j> .  Furthermore,  the  eccentricity  and  the  maximum  possible  electrodynamic 
torque  for  a  given  tether  design  are  both  small,  i.e.  e«l,  s  «  1  .  This  method  would  work  with  eccentric 
orbits  as  well,  but  in  deriving  the  asymptotic  expansion  for  the  libration  angle,  one  would  need  to  expand 
about  the  difference  from  the  reference  eccentricity.  Because  an  EDT  must  orbit  low  enough  to  take 
advantage  of  the  Earth’s  magnetic  field,  the  orbit  is  nearly  circular  by  necessity,  so  the  problem  posed  here 
is  for  a  nearly  circular  orbit.  The  optimal  control  problem  is  constructed  in  a  manner  similar  to  the  ones 
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posed  in  Chapter  III,  but  with  the  additional  constraints  on  the  mean  square  libration  state,  z  ,  and  may  be 
written  as  the  following. 

Minimize  Cost:  J  =  tf 

Subject  to: 

%  =  f(X’U) 

e0  (x(r0))  =  [a0,e0,/0,z0]r  =  [6648  km,0.005,30°,0.0038]r 
ef(x(7}))  =  [a/,//,z/J  =[6658  km,30.5°,0.0014f 
&  («(^))  =  C -2-25  <0  Amps2 

where  x(T)  =  [a,h,k,i,z]T  represents  the  average  states  with  averaged  dynamics 
f  (x(r),u(r))  a  Ax/A T  that  are  described  by  Eq.  (51)  and  Eq.  (11).  Box  constraints  in  Eq.  (17)  are  still 

enforced  as  well  and  the  rms  current  is  defined  by  Eq.  (14).  The  500  kg,  4  km  tether  in  this  case  is  modeled 
using  a  230  kg  upper  end  body,  and  a  220  kg  lower  end  body  with  the  same  current  properties  as  the  tether 
modeled  in  the  Chapter  III  examples.  Solving  this  problem  using  DIDO  for  the  no  drag  case  yields  an 
optimal  control  solution  (Figure  27)  that  drives  the  libration  magnitude  to  the  final  desired  value  while 
executing  the  desired  orbital  maneuver  (Figure  28).  Similar  to  the  minimum  time  solution  obtained  in 
Chapter  III,  much  of  the  thrust  is  used  to  achieve  the  inclination  change  through  the  u4  and  u5  control 
coefficients  corresponding  with  the  frequency  components  twice  that  of  the  orbital  frequency.  In  this  case, 
however,  there  is  a  small  component  of  the  periodic  current  allotted  to  u2  and  m3  to  drive  the  libration 
amplitude  to  the  desired  final  state.  The  libration  angle,  $ ,  depicted  in  Figure  28  was  propagated  using  the 
exact  equations  of  motion  given  by  Eq.  (29)  with  a  stiff  ode  solver  (Matlab’s  ode23t)  to  verify  the  accuracy 
of  the  dynamic  model  and  the  assumptions.  Control  current  is  constantly  phased  with  the  librational 
motion  to  account  for  a  small  frequency  shift  due  to  numerical  errors  in  the  ode  propagation,  i.e. 

k=v(  i+<r) 

where  vc  is  the  true  anomaly  argument  used  in  the  clock  time  domain  controller  given  in  Eq.  (1)  and  C  is  a 
small  parameter  determined  by  observing  the  errors  incurred  during  propagation  of  the  homogenous 
solution  to  Eq.  (46)  (see  Appendix  H).  The  DIDO  solution  shows  the  z  state  history  transformed  into  a 
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Control  Fourier  Coefficients,  u(T) 


Figure  27.  Control  for  Minimum  Time  Orbit  Change  with  Libration  Control,  No  Drag 

libration  angle  history  which  forms  an  envelope  for  the  rapidly  varying  libration  angle.  The  propagated 
orbital  trajectory  and  maximum  libration  angle  envelope  matches  well  with  the  propagated  values 
indicating  that  the  proposed  model  is  sufficient  for  this  problem. 
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Minimum  Time  Orbit  Change,  32  Nodes 


Figure  28.  Minimum  Time  Orbit  Change  State  Trajectory,  No  Drag.  Stars  indicate  DIDO 
derived  libration  envelope;  lines  indicate  propagated  instantaneous  state. 


For  comparison,  a  similar  constant  eccentricity  optimal  maneuver  was  executed  without  any 
restriction  on  the  libration  mean  square  state.  The  constraints  in  Eq.  (52)  were  enforced  with  the  following 
exception  and  addition. 

ef(x(7}))  =  [a/,//]r=[6658  km,  30.5°  f 
g2(u(r))  =  /r+£2-e2=0 

The  additional  path  constraint  is  written  to  enforce  a  constant  eccentricity  maneuver.  The  resulting  control 
profile  (Figure  29)  and  trajectory  (Figure  30)  demonstrate  that  the  maneuver  is  only  marginally  quicker  (by 
a  single  revolution)  but  the  libration  amplitude,  left  uncontrolled,  remains  practically  unchanged  for  this 
time  span.  Given  enough  time,  however,  this  amplitude  can  grow  in  a  thrusting  tether,  so  it  is  important  to 
manage  the  libration  while  maneuvering  an  EDT.  This  is  especially  true  for  a  tether  that  is  long,  carries  a 
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large  control  current,  or  has  a  large  mass  differential  between  the  upper  and  lower  masses  resulting  in  a 
large  electrodynamic  torque  when  the  EDT  is  active. 


Control  Fourier  Coefficients,  u(T) 
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Figure  29.  Control  Profile  for  a  Minimum  Time  Orbit  Change  with  No  Libration  Control,  No 
Drag 

Including  a  state  in  the  dynamics  that  captures  the  magnitude  of  the  out-of-plane  tether  libration 
provides  a  higher  fidelity  constraint  model  that  enables  more  accurate  optimal  control  of  an  EDT.  Since 
the  long  time  scale  equations  of  motion  for  orbit  transfers  assume  a  near  nadir-pointing  tether,  bounding  the 
libration  is  even  more  critical.  The  results  of  this  section  demonstrate  that  it  is  possible  to  control  tether 
libration  while  simultaneously  maneuvering  to  a  new  orbit  using  periodic  control  of  the  EDT  current  over  a 
long  time  scale. 
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Minimum  Time  Orbit  Change,  32  Nodes 


Revolutions 

Figure  30.  Minimum  Time  Orbit  Change  State  Trajectory  without  Libration  Control 

When  drag  is  included  in  the  dynamic  model  (i.e.  D  +  0  in  Eq.(l  1)),  the  controls  initially  boost  the  satellite 
to  take  advantage  of  the  lower  atmospheric  density  at  higher  altitudes  allowing  more  power  to  be  dedicated 
to  cranking  the  inclination  in  the  same  manner  as  the  example  in  Chapter  III.  The  controls  are  shown  in 
Figure  31  with  the  resulting  trajectory  in  Figure  32.  With  drag,  this  maneuver  takes  three  more  days  to 
complete  than  its  no-drag  counterpart,  requiring  a  total  of  270  revolutions. 

Although  we  do  not  demonstrate  definitive  optimality  of  the  control  solution,  compliance  with  one 
transversality  condition  necessary  for  optimality  is  shown.  Because  there  is  no  explicit  time  dependence  in 

the  Lagrangian  of  the  Hamiltonian  of  this  optimal  control  problem  we  have  H  =  0  .  The  Lagrangian  of  the 
Hamiltonian  was  defined  in  Chapter  III  as 

H  =  H  +  Juggl+filx  +  filu 

where  the  Hamiltonian  is  defined  by  H  =  ll f  given  the  Mayer  cost  chosen  in  this  example  and  t. 
represents  the  costates.  Recall  that  the  covector  functions  associated  with  the  path  constraint,  state-variable 
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box  constraints  and  control-variable  constraints  are  represented  by  /j  ,  p  and  fil(  respectively. 
Furthermore,  since  the  problem  posed  here  is  a  minimum  final  time  problem  we  have  H  [tf )  =  —1 ,  so  we 
have  a  condition  that  holds  throughout  the  trajectory,  namely 

H(t)  =  -\  (53) 

DIDO  uses  the  Covector  Mapping  Principle  to  produce  adjoints  and  the  Hamiltonian  as  part  of  the  solution. 
To  check  this  optimality  condition  we  plotted  the  output  Hamiltonian  and  discovered  that  it  indeed  satisfied 
Eq.  (53)  throughout  the  maneuver  within  a  tolerance  of  0.002. 


Control  Fourier  Coefficients 
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Figure  31.  Control  Profile  for  a  Minimum  Time  Orbit  Change  with  Libration  Control  and 
Drag 
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Inclination,  deg  Altitude,  km 
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Figure  32.  Minimum  Time  Orbit  Change  State  Trajectory  with  Drag 
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VI.  Summary,  Conclusions  and  Future  Work 


This  research  demonstrates  how  using  the  method  of  averaging  and  multiple  time  scales  can  be 
used  to  achieve  optimal  controls  for  systems  exhibiting  periodic  behavior,  such  as  maneuvering  low  thrust 
satellites.  Optimal  control  problems  for  maneuvering  electrodynamic  tethers  were  posed  using  averaged 
state  dynamics  and  constraints  and  then  solved  using  a  pseudospectral  optimization  method.  It  was  shown 
that  some  classes  of  complex  optimal  control  problems  that  use  instantaneous  state  dynamics  requiring 
hundreds  or  thousands  of  collocation  node  points  for  accurate  solutions  in  the  clock  time  domain  can  be 
reduced  to  relatively  simple  problems  in  Fourier  space  using  only  tens  of  node  points.  Using  this  method 
of  large  time  scale  optimal  control,  it  is  possible  to  determine  optimal  solutions  of  nearly  periodic  systems 
more  accurately  and  more  quickly  than  optimization  using  the  instantaneous  states.  For  long  term 
maneuvers  spanning  days,  weeks  or  months,  it  may  be  difficult  or  impossible  to  achieve  accurate  solutions 
that  use  instantaneous  state  dynamics  due  to  numerical  round  off  errors.  Using  averaged  state  dynamics, 
however,  small  periodic  behavior  over  each  orbit  is  ignored  enabling  the  optimizer  to  determine  a  trajectory 
for  the  averaged  state,  thus  optimizing  only  the  secular  behavior.  This  greatly  reduces  the  scale  of  the 
problem  for  the  optimizer.  This  method  of  optimal  control  in  Fourier  space  could  assist  engineers  with 
initial  trade  studies  to  determine  design  and  performance  parameters  for  a  tether  or  any  other  low  thrust 
maneuvering  satellite. 

It  was  further  shown  that  this  method  of  large  time  scale  optimal  control  may  be  adapted  to 
accommodate  dynamics  operating  over  multiple  time  scales.  For  the  electrodynamic  tether  controller 
model,  it  was  necessary  to  include  the  effects  of  a  tilted  Earth  magnetic  dipole  which  rotates  once  per  day, 
slower  than  the  orbital  rate  of  the  satellite  system.  None  of  the  controllers  described  in  the  literature 
addressed  a  tilted  Earth  magnetic  dipole  or  an  atmospheric  drag  model  for  electrodynamic  tethers  in  orbit 
transfer,  so  a  model  was  derived  that  included  both.  The  periodic  controller  was  modified  to  include  terms 
resonant  with  the  Earth’s  rotation  and  more  accurate  results  were  achieved  and  verified  against  a  “truth” 
model. 
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To  provide  an  even  higher  fidelity  controller  model,  optimal  libration  control  was  also  examined. 

It  was  shown  that  a  rapidly  changing  state  such  as  libration  may  be  controlled  in  Fourier  Space  by  defining 
a  mean  square  state  and  using  the  averaged  dynamic  equations  of  motion  as  constraints  in  the  optimal 
control  problem.  In  this  manner,  it  was  possible  to  achieve  minimum  time  orbit  transfers  that 
simultaneously  drove  down  libration  amplitude.  Using  this  method  of  large  time  scale  optimal  control  in 
concert  with  instantaneous  state  controllers  operating  in  real  time  could  enable  maneuvering 
electrodynamic  tether  satellites  to  achieve  long  term  transfers  unachievable  using  instantaneous  state 
controllers  alone. 

Optimal  controls  for  low  thrust  satellites  performing  orbital  maneuvers  using  multiple  time  scales 
is  a  wide  open  field  with  plenty  of  areas  to  be  explored.  The  following  is  a  list  of  recommended  follow  on 
research. 

Apply  the  method  of  multiple  time  scale  optimal  control  to  systems  subject  to  a  different  class  of 
dynamic  equations  of  motion.  In  this  research,  optimal  control  problems  were  reduced  in  Fourier  space 
because  we  exploited  what  we  knew  about  the  problem,  namely  that  the  dynamic  system  had  a  fast  time 
periodic  piece  and  a  slowly  varying  secular  piece.  There  are  many  other  systems  that  fall  into  this  category 
that  could  use  this  method.  Additionally,  we  are  not  constrained  to  periodic  systems.  A  basis  in  Fourier 
space  was  chosen  here  because  of  the  periodic  nature  of  the  orbital  mechanics,  however  a  different  problem 
might  be  better  served  in  a  polynomial  space  with  an  orthogonal  polynomial  basis. 

Demonstrate  a  powered  sun-synchronous  orbit.  Non-thrusting  satellites  may  be  placed  in  sun- 
synchronous  orbits  that  take  advantage  of  the  Earth’s  oblateness  in  such  a  way  that  the  orbital  plane 
precesses  once  per  year.  However,  these  orbits  are  typically  constrained  by  altitude,  inclination  and 
eccentricity.  A  continuously  thrusting  system  however  could  potentially  achieve  otherwise  unachievable 
sun-synchronous  orbits.  One  advantage  would  be  that  a  satellite  could  reside  in  a  desired  orbit  while 
maintaining  optimal  solar  panel  orientation  with  respect  to  the  sun  at  all  times. 

Demonstrate  optimal  control  using  a  higher  fidelity  model.  Other  periodic  effects  may  be 
included  into  the  dynamic  model.  A  diurnally  varying  atmosphere,  by  definition,  differs  on  the  day  and 
night  sides.  The  controller,  as  described  by  Eq.  (1),  already  contains  control  coefficients  to  affect 
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perturbations  resonant  with  the  orbital  motion  such  as  diurnally  varying  phenomena.  When  we  introduced 
the  Earth’s  rotating  tilted  magnetic  dipole  into  the  model,  we  had  an  effect  that  is  not  resonant  with  the  first 
two  harmonics  of  the  Eq.  (1)  controller.  New  terms  were  be  added  to  the  controller  to  accommodate  new 
perturbation  effects  (Eq.  (25)).  There  may  be  other  multiple  time  scale  effects  to  consider  in  the  model  that 
operate  at  different  frequencies. 
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Appendix  A:  Derivation  of  Libration  Equations  of  Motion 


In  developing  an  orbital  maneuvering  controller,  it  is  important  to  understand  the  attitude  behavior 
of  the  tether  motion  when  subjected  to  electrodynamic  and  aerodynamic  forces.  Because  of  the  length  of 
the  tether,  gravity  gradient  restoring  torques  can  be  significant.  In  order  to  make  the  equations  as  general 
as  possible  for  2-ball  tether  designs,  few  assumptions  were  made  with  regard  to  the  mass  distribution.  The 
tether  is  modeled  as  two  end-masses  connected  by  a  straight  tether  in  constant  tension  with  a  uniform  mass 
distribution. 

The  conservative  gravitational  force  plays  a  large  role  in  the  tether  attitude  dynamics  and  lends 
itself  well  to  the  development  of  equations  of  motion  using  the  Lagrangian  method.  Constructing  the 
Lagrangian,  we  need  adequate  expressions  for  both  the  kinetic  and  potential  energies.  Using  the 
coordinates  defined  in  Figure  33,  we  can  write  the  endmass  and  tether  velocities  and  thus  the  kinetic 
energy. 


Figure  33.  Rotating  Frame  Coordinates 
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The  inertial  frame  is  centered  at  the  center  of  the  Earth.  The  rotating  frame  is  located  at  a  position 
r  with  respect  to  the  inertial  frame  and  is  centered  at  the  system  COM.  It  consists  of  three  mutually 
orthogonal  unit  vectors;  er  in  the  zenith  direction,  e,  in  the  transverse  direction  and  en  completing  the  triad 
in  the  direction  of  the  angular  momentum  perpendicular  to  the  orbital  plane.  The  vectors  along  the  straight 
tether  extending  from  the  COM  to  mass  1  and  mass  2  are  p,  and  p2  respectively. 


Kinetic  Energy 


The  kinetic  energy  for  the  system  is  derived  by  summing  the  separate  kinetic  energies  for  the  end- 
masses  and  the  integrated  kinetic  energy  across  the  length  of  the  tether.  The  velocity  for  mass  1  and  its 
square  are 

Vi=r  +  Pi=r  +  p;+(o0+©e)xp1=r  +  p;+£lxp1 

Vj-Vj  =(f+p;  +  QxPl)-(f+p;+QxPl) 

=  r -f  +  2f -pj  +2r -(flxpj  +  pj  -pj  +2p[  •(J2xp1)  +  (flxp|)-(flxp1) 
where  ioe  is  the  angular  velocity  of  the  straight  tether  in  the  rotating  frame,  u>0  is  the  angular  velocity  of 
the  rotating  frame  with  respect  to  the  inertial  frame  and  =  coo  +  ioe .  Primed  vectors  indicate  radial 
derivatives  with  respect  to  the  rotating  frame  (e-frame).  Given  the  length  of  the  straight  tether,  L,  we  may 
write  the  relative  position  vectors  p,  and  p2  in  terms  of  a  reduced  mass  and  L. 
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where  p  is  the  unit  vector  in  the  direction  along  the  straight  tether  from  mi  to  mi.  The  reduced  mass  /im 
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mm 

where  Mx  =  mx  +  —^~,  M2  =  m2+^~  and  M  =  w,  +  m2  +  mt ,  the  total  mass.  Derivation  for  this  3 -body 

reduced  mass  is  found  in  Appendix  B.  Considering  the  appropriate  substitutions  for  p,  and  p2 ,  the  kinetic 
energy  of  mass  1  is 

T  1 

Tl  =-'»iV,  -V, 

2  2M,  [  M,  '  V  J  \  JJj 


where  L'  is  the  time  derivative  of  L  w.r.t.  the  e-frame.  Likewise  the  kinetic  energy  for  mass  2  is 


Note  that  the  fourth  term  in  braces  vanishes  since  L'  •  £1  x  L  =  0  . 

The  kinetic  energy  of  the  tether,  however,  must  be  integrated  from  tip  to  tip  (i.e.  —px  to  p2  as 
shown  in  Figure  34) 

1 

T,  =  ^  f  •  \tdm 

-pi 


Figure  34.  Straight  Tether  Integration 
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where  the  velocity  of  any  given  element  ds  along  the  tether  is 

vt  =  f  +  s  =  r  +  s'  +  flxs 

and  a  mass  element  for  a  tether  of  uniform  density  is 


dm  -  —ds 


Substituting  in  appropriate  terms,  the  tether  kinetic  energy  may  be  written 

1  ^  ft 

Tt  =—  —  j  r-r +2r -s'  +  2r  •Qxs  +  s'-s'  +  2s'-f2xs  +  (i2xs)-(f2xs)tfc 


Recognizing  that  a  tether  section  that  spans  length  s  stretches  in  linear  proportion  to  the  overall  tether 

/  sL  /v 

stretch,  i.e.  s'  =  y^  L'  =  - — p  and  we  can  write  the  kinetic  energy  integrand  in  terms  of  the  scalar  s. 


yyi  i  SJ-j  |*  J_j  /•  |> 

Tt  = — L{f-fZ  +  2r-p  f  - — fifr  +  2f-flxp  1  sds  +  —  f  s2 ds  +  2pflxp  —  1  s2ds 

2L  J  /  j  i  j  j  j 


sL 

L 


riY  « 

L 


XL  , 

\  J 


(54) 


H2 

-(flxp)-(Qxp)  j  s2ds} 


u  u. 

Since  -/?,  = - —L  and  p2  =  — —L  ,  we  can  rewrite  the  integration  limits  and  the  tether  kinetic  energy 


M, 


becomes 


nn  ,■  •  „  •  ~  L  L 

T.  =  — -  {r  •  rL  +  2r  •  p - 

2  L  12 


2  (  ..  2  ..  2  j2 


Mm  M, 


km2  m2  j 


+  2r  Qxp — 
2 


2L2\ 

(  3  3  \ 

Mm  +  Mm 

3  1 

M  3  M3 

V 1VL  2  1V1\  J 

-(£lxp)-(£lxp)y 


(  M,n  ,  M, 


V^V  M,  , 


(  Mr,"  ^ 

yM2  M,  , 

3 

1} 


Where  the  fifth  term  in  the  braces  of  Eq.  (54)  drops  out  since  p  •  Q  x  p  =  0  . 
Thus  we  are  left  with 


nn  ^  .  „ 

T,=  —rr  +  m,Lrp 


2  M 


-m?Zr-12xp  — - 


jj-  f  a  j  3 


6 


M{  +M2 


M 


3  'A  m  r 2  (  ,^-3 

(flxp)-(flxp)— 


2  M  J 

3 


M{  +M2 


M 


J 


Note  that 


80 


2 


2 


fim2  M,,,2  _  -  M2  _  m\  +  mimt  -  m2  ~  m2mi  _  imi  -  m2 )  ( m\  +tn2+  m, )  _  ml  -  m2 

Mj-JwJ-  W  ~~  M2  ~  M2  M 


The  total  kinetic  energy  of  the  system  is  then 


T  =  Tl+T2+T, 


1  , .  ~m,  m2  1  m,  ( m,-m2 

=  -Mr-r  +  <um  L  +  ^-  +  —  — - - 

2  M2)  2  y  M 


(Zr-p  +  Tr  Qxpj- 


+m,>  1  U  ,  (nxfi).(nx^ 

2 2 6  I  M1  I  '  '  'y  " 


The  term  in  the  first  set  of  braces  vanishes  as  shown  below. 


M l  M,  +  M2J  2  t  M 


-M 2ni\  +Mlm1  +  -  m2 ) 


(  m, }  (  \  m,  t  \ 

”i 1  m2  +—\  +  m2 1  m<  +—\  +  —'y m'  - m2 ) 


Furthermore,  the  quantity  (fix  p)  •  (fix  p)  is  actually  Q  JQ  or  fi  •  J  •  fi  in  matrix  and  tensor  form 


respectively  where  J  is  the  specific  inertia  tensor  (see  Appendix  C).  Thus  the  total  kinetic  energy  for  the 
system  is 


„  1,,.  ■  if  ->(  m\  w,  1  m,(  M.3  +M23]\  ( -2  ,  -  ^ 

T  =  — Mr -r +  — < /j  2  +  — 1 — -2-  \\l}  +LSI- J  il 

2  2  \m  {M2  M2)  3  ^  M3  Jj\ 


The  mass  term  in  braces  may  be  shown  to  reduce  to  an  equivalent  reduced  mass  fie  as  shown  in  Appendix 
B. 


f  \ 

ml  m2 

m, 

_i _ L 

' m3+m2 n 

mt 

[m2+  M2) 

i 

3 

M3 

v  1V1  / 

~  f^m  ~7~ 
O 

The  final  form  of  the  total  kinetic  energy  is 


T  =  -Mr-r  +  -u  (z2+Z2njft 
2  2  e\ 


81 


The  first  term  on  the  right  hand  side  of  Eq.  (55)  is  the  translational  energy  of  the  whole  system  acting 
through  the  system  COM.  The  second  term  accounts  for  the  rotational  energy  acting  through  the  COM  and 
the  relative  motion  energy  between  the  two  end-masses. 

Potential  Energy 


Figure  35.  Position  of  Endmass  1 

The  potential  energy  is  derived  for  each  end-mass  at  its  distance  from  the  center  of  the  Earth.  The 
potential  energy  of  the  tether  wire  is  integrated  for  each  elemental  mass  along  the  length  of  the  wire.  To 
obtain  the  potential  relative  to  the  COM,  a  binomial  expansion  of  the  potential  energy  expression  is  used. 
For  mass  1,  the  potential  energy  may  be  expressed  in  terms  of  the  radius  vector  to  the  COM  of  the  system, 
r  ,  and  the  n^  position  vector  relative  to  the  COM 

-ngmx  -/y«. 


V,=- 


1  +  2^+<\ 


2 


where  we  have  used  the  geometry  in  Figure  35  and  the  law  of  cosines  to  infer  that 


r,2  =  r2  +  pt 2  -  2rpl  cos  6  =  r2 


'i+4+2fcV 


\  P\  2Pr02 

2 


An  expression  for  rx  1  is  derived  using  a  binomial  expansion  and  ignoring  terms  greater  than  or  equal  to 


0 


P~ 
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r  2  r  2  r 


Thus  the  potential  energy  of  mass  1  is 
vi  =  - urn ,/f1 


^l<;nh  L  P,-er  3  f  p]  er  V 

r  I  r  2 1  r  J  2  r2 


Substituting  in  the  expressions  for  p,  and  p2  we  obtain 


~Msr>h 

1^  x  V,nL  1 

f  P  '  O 

,3^1 

(  P'^ 

I2  1/C£21 

r 

l  «. 

l  r  J 

2  M;  1 

l  r  X 

1  2  M{r 2 

^Wijl  +  ^(p-er)  +  ^S 

1  M,r  7  2  M{r2 


(3(P-er)2-l)} 


Similarly,  the  potential  energy  of  mass  2  is 


F, 


M2rV  7  2M,V 


(3(per)2 


-1 


The  potential  energy  of  the  tether  is  the  integration  of  all  the  elemental  potential  energies  along  the  tether 
length. 


y, =-»,! 


- ^—^-d m  where  dm  =  —ds 

r  +  ser  L 


Using  a  binomial  expansion. 
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The  total  potential  energy  is  then 
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V  =  Vl  +  V2  +  V, 

„  VSM  MgL  nmm2_+m1(m1-mx 

2 


r 

ST 

I'Y 

1 

^ \ 

|  M22m1  +  M2m-, 

2  r3 

M2 

(  \/t  3 


3 


Mj  +  M. 
M3 


(P'®r)‘ 


(3  (p  •  er )“  -l) 


The  first  term  in  braces  vanishes  to  zero  as  shown. 


The  second  term  in  braces  is  the  equivalent  reduced  mass,  fie .  The  total  potential  energy  is  then 


2r 3 


^(3(p-er)2-l) 


(56) 


The  first  term  in  (56)  is  the  potential  energy  of  the  entire  system  mass  acting  as  a  point  mass  at  the  COM. 
The  second  term  is  the  gravity  gradient  potential  due  to  the  center  of  gravity  offset  from  the  COM. 


The  Lagrangian  Equations  of  Motion 


With  the  kinetic  and  potential  energies  in  hand  (equations  (55)  and  (56)),  we  may  now  construct  the 
Lagrangian  function,  £,  =  T—V,  and  form  the  Lagrangian  equations  of  motion. 


1 


1 


L  =  -Mv2  +-u,  C2+Z2n-J-il + 


ft/'/ 


2r 


(3(p-er)2-l) 


To  get  the  equations  of  motion  in  terms  of  the  in-plane  and  out-of-plane  libration  angles  of  the  straight 
tether,  we  need  reference  frames  and  coordinates  with  which  to  evaluate  the  vector  operations  to  obtain 
scalars.  The  body  frame  and  the  rotating  frame  (e-frame)  will  serve  well  for  this  purpose.  The  rotational 

energy  term,  L2£l  ■  J  ■  Q  ,  is  evaluated  using  the  body  frame  depicted  in  Figure  36.  The  orbital  rotation  may 
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be  expressed  in  body  coordinates  as  wa  = 


cosix\<j) 

0 

®COS(2$ 


and  similarly  the  tether  rotation  in  the  rotating  frame 


may  be  expressed  in  body  coordinates  as  coe  = 


9  sin 
-<!> 

9  cos^ 


Letting  s9  =  sin  9  and  c9  =  cos  9  we  may  write  the  total  inertial  angular  rate  of  the  tether  system  in  body 
coordinates  as 


fi  =  to0  +<ne  = 


(<9  +  (o)s(/) 


The  specific  inertia  tensor  may  be  expressed  in  body  coordinates  as  well  in  matrix  form.  Approximating 
the  tether  as  a  thin  rod  with  zero  inertia  about  the  b,  axis,  we  write 


SI;  = 


0  0  0 


0  L- 
0  0 


0 


thus 


L-n  J  n  =  Z/Q'  JQ  =  L-  w  +  c-<t>  9  +  CO 


^ 6 2  +c2<f>[9  +  co}  j 
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where  co  is  not  constant  in  general  and  is  equal  to  the  rate  of  change  of  the  true  anomaly,  v  .  To 
evaluate  (  p  •  er )  ,  the  e-frame  is  convenient.  In  the  rotating  e-frame,  the  length  vector  and  its  derivative 


may  be  written 


L  = 


L 


cos  ^  cos  9 
cos  <!>  sin  9 
sin^S 


=  Lp 


cos  ^  cos  9 

-9  cos  ^  sin  9  -  $  sin  ^  cos  9 

L  =  L'  +  ©e  x  L  =  ip  +  Zp  =  L 

cos  (/>  sin  9 

+  L 

9  cos  (/>  cos  9  -  (j>  sin  (/>  sin  9 

sin  (f) 

e 

</)COS(/> 

c<j>c9 

c(/>s6 

S(j) 


,  p  ■  er  =  c(/>c9,  (p  •  er  f  =  c2(/>c29 


The  Lagrangian  is  then 

£  =  Imv2  +  |Z2  +  i2  +  cV(6*  +  v>)2 ) 


c~0\d  +  vS\\  +  -^—+  [3c2(j)c2e-\) 


2  r 


The  in-plane  libration  equation  of  motion  may  be  obtained  as  follows. 


dIi  =  jueL2c2</>(0  +  v),  —  f^]  =  iueT2cV((9'  +  u)  +  /«e((-2Z2c^#)  +  2TZc2^)(i9  +  v) 


d6 

dL  -3jUeMeL2 

dd 


dt  v  d6 

gt~e  c2  <j>cds6 


dt 


dL  'j  dL  _ 

^de)~de~Qe 


( 

c2^(#  +  i>‘)  +  2 

V 

jc2<!>  -  0c<j>s<])\j^9  +  v^  +  -^c2<j>c9s9 

ii 

The  out-of-plane  libration  equation  of  motion  is  similarly  obtained. 

=  -jueL 2  c</>s</>(e  +  vy  +-^-c<j>s<j>c29 

dL  '  dL')  dL  _ 
dt  Kd(j>)  d</>  * 
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(57) 
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(58) 
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l  L 

('M‘  +  7 r‘le 

V  A 

\ 

/ 

The  motion  along  the  length  of  the  tether  is  free  when  the  tether  is  slack,  but  constrained  when  the  tether  is 
in  tension.  This  equation  of  motion  is  given  by 


dt 


dL  I  T 
— -  =  /f.T 

8L  1  e 


—  =  MeL' 

8L  e 


dt 


dL )  dL 


8L'  J  dL  Ql 


(jr  +  c2<f>{e  +  v'j  +^(3c2</>c20-l) 


VeL 


--if  -cV((9  +  t>)2  (3cVc26>-l) 


=  Ql 


(59) 


To  complete  the  right  sides  of  these  equations  of  motion,  we  need  to  write  the  non-conservative  generalized 
forces  acting  on  the  system. 


Non-Conservative  Generalized  Forces 

There  are  two  external  non-conservative  forces  considered  in  these  equations  of  motion, 
aerodynamic  drag  and  Lorenz  force  due  to  the  current  in  the  wire  interacting  with  the  Earth’s  magnetic 
field.  At  the  operational  altitudes  of  interest,  the  magnitudes  of  both  of  these  forces  cannot  be  neglected. 

In  this  section  we  will  derive  the  generalized  forces  on  the  tether  due  to  these  effects. 

Aerodynamic  Drag 

Although  the  air  density  is  very  low  in  the  stratosphere,  drag  is  not  negligible.  Over  many  orbits, 
the  atmospheric  drag  will  eventually  cause  the  tether  orbit  to  decay  if  there  is  no  restoring  force.  The  air  is 
too  thin  to  model  as  a  fluid,  i.e.  the  molecular  mean  free  path  is  large  compared  to  the  dimensions  of  the 
satellite.  Therefore  we  use  a  free-molecular  flow  model  and  only  consider  a  drag  force  opposite  the 
direction  of  the  velocity.  This  force  acts  on  both  end-masses  and  the  tether  itself.  Each  end-mass  has  a 
different  ballistic  coefficient  and  the  system  COM  is  not,  in  general,  in  the  center  of  the  tether. 
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Furthermore  the  atmospheric  density  varies  exponentially  along  the  length  of  the  tether,  thus  the  impact 
force  of  incoming  atmospheric  particles  varies  with  height  along  the  wire  as  shown  in  Figure  37.  Due  to 
the  uneven  distribution  of  aerodynamic  forces  about  the  COM,  there  will  be  aerodynamic  torque  acting  on 
the  tether  system. 


but 


Figure  37.  Tether  Subject  to  Atmospheric  Drag 


The  generalized  in-plane  aerodynamic  torque  is  given  by 


_  5vi  ^  Sv2  ..  dfl 

Q„  —  F, —  +  F, —  +  M.  — — 

6a  de  do  dd 


dv,  ^(f  +  Pj)  dp,  5(p;  +HxPl)  dpi  ~  0p.  dfl  dfl 

— F  =  — — ; — -  =  — F  =  — - ; - -  =  — f  +  llx — F  + — -xp.  = — -xp. 
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since  — r  =  — F  =  — F  =  0  .  Recalling  that  Q  =  to  +  to  ,  we  have 

de  de  de  " 


5vj  3(to  +to  )  0oe  .  0too 
— F  = - ; - -xp  = — ^xp  since  — F-  =  0. 
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Similarly 


— I-  =  — ^  xp,.  thus  the  aerodynamic  torque  reduces  to 
86  86 


Qea 


xp,  +F2  • 


dtoe 

66 


xPi  +Mt  • 


oil 


The  forces  on  the  end-masses  depend  on  the  respective  ballistic  coefficients  and  air  densities  and  vehicle 
velocity.  For  the  purpose  of  computing  the  aerodynamic  forces  on  both  end-masses,  the  velocity  due  to 
orbital  motion  is  considered  far  greater  than  the  velocity  of  each  end-mass  relative  to  the  COM  due  to 
rotation.  Also,  the  tether  is  considerably  shorter  than  the  distance  to  the  center  of  the  Earth,  so  the 
directions  of  F,  and  F2  are  approximately  parallel  with  the  velocity  vector  of  the  COM.  In  the  force 


computations,  we  use  Va  as  the  velocity  of  air  molecules  impinging  on  the  vehicle.  Ignoring  the  rotation 

of  the  atmosphere  since  this  velocity  contribution  is  miniscule  compared  to  the  orbital  velocity  in  low  Earth 
orbit,  we  can  write  the  air  velocity  in  the  rotating  e-frame  (designated  with  a  subscript  “e”)  in  terms  of  the 
flight  path  angle  y  (Figure  38)  as 


-sin/ 

where  va  = 

-sin/ 
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Thus,  the  aerodynamic  force  on  mass  1  is  described  as 

-sy 
-cy 
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CdiAPoe 


~(*+Pi  er 
A* 
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~  CdiA\ pae 


where  h  is  the  altitude  of  the  COM,  h  is  the  characteristic  height  of  the  atmosphere,  wp  is  the  argument 
of  perigee,  Crfl  is  the  effective  coefficient  of  drag  on  M] ,  At  is  the  presented  area  of  mass  1,  and  /90isthe 
atmospheric  density  at  the  COM  altitude.  Also 


de  M 


-af)s6 

do  /umL 

and  — ^xp,  =  —H— 
86  M2 
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The  torque  about  the  perpendicular  to  the  orbital  plane  about  the  COM  due  to  the  drag  on  mass  1  is 
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( -sysd  +  cycQ ) 
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=  ~BiPoe  h 
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where  a*  = 

Mi 


is  the  effective  ballistic  coefficient  with  Mx  in  the  denominator,  not  mx .  Likewise,  the 


aerodynamic  force  on  mass  2  is 
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with  a  corresponding  aerodynamic  torque  of 


-f  h+^a/icg) 
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The  moment  due  to  the  varying  air  drag  on  the  tether  must  be  integrated  along  the  tether  length.  An 
infinitesimal  force  acting  on  an  element  of  surface  area  on  the  tether  is  proportional  to  the  presented  area  to 
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the  air  flow,  the  velocity  squared,  the  coefficient  of  drag,  and  the  air  density  at  that  position  on  the  tether. 

Its  direction  is  in  the  direction  of  the  air  flow. 

dF  =C^p(h)va2dAvli 

The  elemental  area  presented  to  the  air  stream  is  dA  =  cos  adtds  where  a  is  the  angle  of  attack  as  shown 
in  Figure  38  and  dt  is  the  tether  diameter  (tether  wire  modeled  as  a  long  cylinder).  Using  a  coefficient  of 
drag  for  an  infinitely  long  cylinder  of  approximately  two  (i.e.  Cit  ~  2  ),  we  may  express  the  moment  on 
tether  alone  due  to  drag  as  (see  Ref  51  pp  250-251) 

ft  ft  ft 

M,  =<|j)sxiiF  =  j  (p  x  va  )sdF  -  (p  x  va )  j  sp{hs)v2dA  =  (px  va)v02  cosadt  J  sp(hs)ds 

-a  -ft  -ft 

where  hs  is  the  altitude  of  the  element  which  may  be  written  hs  =  (/?er +.sp)-er  =  h+scost/>cos0 .  Fora 

more  thorough  discussion  of  the  tether  coefficient  of  drag,  see  Ref  52.  Integrating  in  this  manner  presumes 
that  cos  a  >  0  since  a  negative  value  would  imply  force  acting  on  the  backside  of  the  elemental  area,  thus 
only  the  area  presented  to  the  air  stream  is  considered. 


Figure  38.  Tether  Element  and  Drag  Geometry 
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The  angle  of  attack  is  defined  as  the  angle  between  the  unit  inward  normal  to  the  elemental  area,  n ,  and 
the  air  stream  velocity  unit  vector,  va  as  shown  in  Figure  38.  Using  the  angles  defined  in  this  figure  we 
can  determine  the  angle  of  attack  in  terms  of  the  libration  angles,  9  and  </; ,  and  the  flight  path  angle,  y , 
using  the  following  spherical  trigonometry  relation. 

cos  P  =  cos  (f>  cos  8  +  sin  <j>  sine?  cos  =  cos^cos<? 

r) 

cos  p  =  cos  (/>  sin  ( 9  +  y) 

n  A  A  A 

P  =  —  —  a  since  n,  va  and  p  all  lie  in  the  same  plane  and  p  _L  n 

.'.  sina  =  cos/?  =  cos^sin(6l  +  ^) 

.  .1/  .  .1/ 

cosa  =  (l-sin2  aj1  =  (l-cos2  (^sin2  (9  +  y))'2 


s=--e-y 

2 


cost? 


{  n 

=  cos - 


u 


(9  + 


Notice  that  for  a  circular  orbit,  y  =  0  for  all  time.  Additionally,  if  there  is  no  out  of  plane  motion 
then  cos  or  =  cost? . 

Air  density  is  modeled  as  an  exponentially  decaying  atmosphere  with  characteristic  height  (or 
scale  height)  h  .  Letting  p0  be  the  reference  density  at  the  altitude  of  the  center  of  mass,  then  the  density 
at  the  tether  element  location  is 


-(/t+SCOS^COS#)  —S  COS  (j)  cos  0  —s  COS  <j)  cos  0 

p{hs )  =  pQe  *’  =  pae  'h’e  h’  =  p(h)e  h’ 

-y, 

where  /?(/?)  =  pae 

Substituting  these  values  into  the  tether  moment  equation  yields 

_hy  |  .  COS  9  COS  j 

Mt  =  v2 p0e  'h  dt  (l-cos2  ^sin2  (&  +  y))  2  (px  v,)  J  se  h  ds 


Mt  =  v2p[h)dt  (l-cos2  f^sin2  (9  +  y)^)  2  (px  va) 


(  ti  V  =^f  C(j>cd  , 
e  " - s-1 


v 


c(j>c6  J 


-cdc<f>pmL 


=  v2  p(ti)c\e 


h  M-) 


-cOc</>pmL 

v  h*  M2 


-1 


\  +c0cfnmL 

h*M, 

-e  1 

c9ct/pmL  ,)} 

/ 

l  h  M,  Jj 

'(Pxv.) 


where  C  =  d,  (l-cos2  (^sin2  [9+ 


c(/>c9 
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Now  the  torque  effect  on  6^  due  to  the  aerodynamic  drag  on  the  tether  alone  is 


d£l 

— --M  = 
dd 


-ct)c(pfimL  f  \  f  \ 

PWcle  ^7^-1  -e  *■“'  1  (d(^9-cyc«)) 


0 

(/il 

where  12  =  —d>cd  — -  =  0  and 

dd 

v  +  d  1 


Mt  =  v2/?(/?)C l  e 


-c0cfjtfimL  f  „  ,  \  +c9cj>fimL  f  n  , 

h’M2  -cOc<j>nmL  1  g  iX  cQc<piimL 
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-sys<f> 
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er 

e, 

K 

cys<j> 

c<j)cd 

c<j>sd 

S(j) 

= 

-sysef) 

-sy 

-cy 

0 

efi^sysd  -eyed) 

The  total  aerodynamic  drag  torque  affecting  #  is  therefore 


^  ^  56,  _  56  d£i 

On  —  F, - —  +  Ft - —  +  M, - — 

dd  dd  dd 


Qea=v2p(h)c<f>(cycd-sysd)l,^-  Bie  M,h  -  B.,e 


Hm  Loped  -  [umLc<fc6 

M,h*  n*  MS 


^{-cdc^mL 


-1  -e 


^(cdc^L 


Notice  that  this  expression  reduces  significantly  if  the  assumption  that  the  atmospheric  density  does  not 
vary  significantly  across  the  length  of  the  tether.  Looking  at  the  expansion  of  the  exponential  terms  in  Eq. 

(60)  and  letting  px  =  an(j  p  -  we  wrjte  for  a  cjrcu|ar  orbit  (  y  =  0  ) 

Myh  "  M2h 


Qea  =  v2p(h)c(f)cd< B‘  1  +  x+y  ~K  l-P2+^~  ~ 


d,  (l_c2  </>s2  dy  l- Pi  +Jy  l  +  Pl+-^~  (p,  - 1) 

\C(pCU )  l  l  l  )  \  L 
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=  v2p ( h)c</>cO  ^y~{bI  (1  +  p1 ) -  B\  (l - p2 )) - d,  (l -c2  0  s2  O)' 


Qea  =  v2p(h)c</>c6 |^y1(-si  (1  +  P\ )~ K 0 ~ Pi)) ~ 
d,(  1-C2,zis2  e)' 


f  ti  ' 

a/>cO 


(. p\  -  A2) 


(  h*  ^ 

2 

jYqzSc#V  2  2 

fill 

yc<j>cO  j 

(l  *•  J '  ■ 

1  ^  — 

The  last  term  in  the  braces  is  rewritten  as 


p2Is 

r*m 


'  1  1  ' 

'  ( m,  -  tn2 )  [} 2 

M\  M2  j 

{  M  ) 

=  "2  P,L 


where  we  have  used  the  mass  term  equalities  given  in  the  kinetic  energy  discussion  in  Appendix  A  and  the 
definition  of  the  COM  offset  distance  p:  given  in  (66).  Therefore  the  in-plane  torque  is 

Qea  =  p{h)v2Lc<j>ce\^-{B\(\  +  px)~ 5*(l-  p1))  +  2dlpl(i-c2  ^s2  tf)7 


which  is  accurate  for  J  «  1  •  This  assumption  would  mean  that  the  length  of  the  tether  is  much 

smaller  than  the  characteristic  height  of  the  atmosphere.  For  altitudes  up  to  about  160  km,  h  «  1km  . 
Between  160  and  400  km,  h  ranges  from  30  to  60  km.  If  we  desire  a  first  order  approximation,  namely 


that  —  <K  1 ,  then  the  torque  due  to  the  tether  drag  would  be  derived  as  follows. 
h 


Qga  =  v2p(h)c<j>cO \^{b; -Bl)-d,  (l-c2  ^s2  of 
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\ 


Pm  ,  P, 

yM2  Mx  j 


recognizing  that  pn 


V  ^2  J 


=  1. 


And  the  total  aerodynamic  drag  torque  with  this  simplifying  assumption  is 


Qea=v2p(h)Lc^cO\^{B; -Bf-^^-c2  ^  0)} 


(62) 
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To  derive  the  torque  affecting  the  out-of-plane  libration  angle  </;  we  follow  the  same  process. 


^  _  3v,  _  5v,  d£2 

Q ,  —  F.  •  — —  + 1'  .  •  — : — h  M  •  — — 

*  d(j>  -  d(j)  d(j> 


,  0v,  56.  uL 

where 

d</  d</>  M, 

The  moments  are 
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The  total  moment  about  the  COM  affecting  (j)  due  to  aerodynamic  drag  is 


Qta  =v2  p{h)s(j)(sycO  +  cysO)\ 

f  -c6c(fnmL  / 
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When  the  tether  is  significantly  shorter  than  the  characteristic  height  of  the  atmosphere  (  2r  <3C  1  | ,  the  air 

v  h 


density  is  considered  constant  along  the  tether  and  the  following  simplifications  can  be  made  using  the 
same  procedure  that  produced  Eq.  (62). 

Q(pa  =v2p(h) s</>(syc9  +  cysG){ ( B\  -B\)-C ( 


1  1 

yM2  Mx  J 
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£?*,  =  v2p(h)s</(syc0  +  cys0)^-^-(B*2  -B{)-dt  (l 
Furthermore,  for  circular  orbits  y  =  0  for  all  time  and  we  have 


Qta  = v  V  ih )  ^  (b;  -  b;  )  -  d,  (1  -  c 


(63) 


As  we  shall  see  later,  this  torque  affecting  out-of-plane  libration  vanishes  near  the  equilibrium  point  where 

e  ->  o,  <t>  ->  o . 

The  generalized  force  due  to  air  drag  along  the  tether  length  is  a  force  given  by 


dL  dL  8L 


d£l 

Since  — —  =  0  the  last  term  drops  out.  The  remaining  two  are  determined  by 
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When  W  <g  1 5  then  the  equation  reduces. 
/  n 


QLa  =^-v2p(h)c0(syc0  +  cys0)(B *  -5*) 

If  the  tether  is  in  a  circular  orbit  such  that  y  =  0  then  the  equation  further  reduces  to 

Qu=!Y*lp{h)ctse(B'l-B'2) 
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Lorenz  Force 


When  there  is  a  current  driven  through  the  conducting  tether  wire  in  a  magnetic  field  a  force  is 
created  in  the  direction  mutually  perpendicular  to  the  local  magnetic  field  lines  and  the  current  direction 
according  to  the  relation  d¥e  =  IdL  x  B  where  /  is  the  current,  B  is  the  Earths  magnetic  flux  density  and 
L  is  the  length  vector  along  the  direction  of  the  straight  tether.  For  a  uniform  current  in  a  straight  line 
electromagnetic  tether  we  simply  have  Fe  =  7L  x  B  .  Using  a  non-tilted  dipole  model  of  the  Earth’s 
magnetic  field,  we  can  write  the  magnetic  flux  density  direction  as  a  function  of  true  anomaly,  v  ,  magnetic 
inclination,  i ,  and  the  argument  of  perigee,  cop  .  The  magnitude  depends  on  the  magnetic  dipole  of  the 

Earth,  Ym  I  *  - Atm3  ) ,  and  the  radial  distance  from  the  Earth’s  center,  r  .  The  magnetic  flux 

V  A  ■  m  ) 

density  is  modeled  here  as 


-2  sin 

+  v)sin/ 
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On  the  surface  of  Earth,  B  =  = 


8e6 
63  783 


=  3e”5 
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■  m 


.  Another  250  km  up  from  the  surface,  the  value  is 


about  2.7 e  5 


N 
A- m 
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A  simplifying  approximation  is  that  the  length  of  the  tether  is  significantly  shorter  than  the  distance 
between  the  center  of  mass  and  the  center  of  the  Earth  <s  1  j .  This  means  that  the  magnetic  field 

strength  and  direction  is  approximately  the  same  across  the  length  of  the  tether.  The  model  is  then  simply  a 
uniform  force  distribution  along  the  tether  as  shown  in  Figure  39.  For  analyses  in  which  we  have  already 
assumed  the  length  to  be  short  compared  to  the  atmospheric  scale  height  (40-50  km)  when  modeling  the 
atmospheric  drag,  this  approximation  is  justified.  The  resultant  force,  Fe ,  may  be  used  to  determine  the 

( m2  —  m, )  L 

moment  about  the  COM  with  moment  arm  p,  =  - - - — p  .  The  force  is 

2  M 

Fe  =  7L  x  B  =  ILp  x  B 
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Bn 

B'Cficd  -  Brc(/>s6 

The  torque  affecting  the  6  state  is 


ftt=F.-T7  =  F' 
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The  torque  affecting  </;  is 
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(65) 


There  is  no  tension  in  the  wire  created  by  the  Lorenz  force  since  the  force  is  perpendicular  to  the  wire  so  no 
work  is  being  done  to  move  the  end-bodies  toward  or  away  from  each  other. 
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aj)c9 
ccj>sO 

s<t>  i 

Qu  =P'Fe  =0  since  plFe 

The  total  generalized  forces  (torques)  are  the  sum  of  those  due  to  the  electrodynamic  and  aerodynamic 
forces. 


5v  d\  dlp'  +  flxp)  dp'  d(Zp) 

Q,  =  —  •  Fp  where  —  =  — - : - -  =  — -  = - =  p  = 

Le  dL  ‘  dL  dL  dL  dL 


Qg  Qga  +  Qge 

Q<j)  Q-t])a 

Ql  =  Qlu  +  Qlr 

Now  these  values  may  be  substituted  into  the  right  hand  sides  of  the  equations  of  motion  (Eqs.  (57),  (58) 
and  (59)).  For  purposes  of  this  research,  the  electrodynamic  tether  will  be  controlled  using  the  current  in 
the  wire. 
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Appendix  B:  Reduced  Mass  Derivation 


In  order  to  capture  the  mass  of  the  tether  and  end-masses  of  the  “dumbbell”  tether  model,  a  reduced  mass  is 
used  to  simplify  the  equations  of  motion.  Starting  with  the  geometry  depicted  in  Figure  40,  we  can  write 
the  positions  of  px ,  p2  and  pt  relative  to  the  COM. 


Figure  40.  Tether  Mass  Distribution  Geometry 


Calculating  the  distances  from  each  end-mass  to  the  COM  we  have 


T 

m, — hm9L 
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U  AA  mi  A  A  mt  AA  .  MM 

where  =  mt  +—,  M2  =  m2  +  — ,  M  =  ml  +  m2  +  m,  and  pm  =  — — — . 

The  distance  from  the  system  COM  and  the  tether  COM  located  at  the  midpoint  of  the  tether  is  given  by 

L\  (V 
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m. - +  m2 , 

11  2  \2 
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2  M 


(66) 


The  equivalent  reduced  mass,  pe ,  given  in  equations  (55)  and  (56)  may  be  reduced  as  follows. 
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Notice  that  if  we  approximate  the  tether  mass  as  zero  we  are  left  with  the  more  familiar  two  body  reduced 
mass. 


i“.  = 


TOjTO2 
TO,  +  TO2 


Many  analysts  prefer  to  model  the  system  as  a  dominant  end-mass  with  all  other  masses  negligible.  An 
example  is  the  space  shuttle  or  space  station  with  m, ,  m,  <K  to,  .  In  this  case,  mass  terms  become 


,,  ,,  m,m, 

=  —77^ 
M, 


.  ,  mt 

=  M,  =  TO,  H - - 
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and  the  equivalent  reduced  mass  reduces  to 


TO,  TO,  TO, 

«  =  TO,  + - =  TO,  - 1 - 

2  6  '3 


102 


Notice  that  m1  H — 


m, 


L2  is  the  moment  of  inertia  of  a  mass  on  the  end  of  a  rigid  rod  about  an  axis  through 


the  other  end  of  the  rod. 

It  may  be  desirable  to  simulate  the  tether  using  finite  elements.  The  elements  may  be  modeled  as 

straight  bar  links  with  no  end-masses,  only  tether  mass.  In  this  case  the  reduced  masses  for  each  element 

may  be  approximated  in  the  following  manner. 

m1  =  m2  =0 

m.  m. 


This  effective  mass,  when  multiplied  by  1 '}  ,  is  the  moment  of  inertia  about  the  center  of  a  uniform  density 
rod. 
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Appendix  C:  Tether  Moment  Of  Inertia 


The  derivation  of  the  specific  inertia  dyadic  is  derived  in  several  dynamics  texts  (see  [53]  for  example),  but 
will  be  repeated  here  for  convenience.  From  a  vector  product  identity,  we  may  write 

=  £n-(p-p)i-n-(fi.p)(n.p) 

=  ^Q-(/72l)-Q-(Q-p)(p-i2) 

=  £n-(/rl)-il-fi-(pp)-n 

=  £1-J-£1 

The  specific  inertia  tensor  is  given  by  J  .  This  tensor  expression  in  matrix  form  is 

(fixp)-(fixp)  =  Qr  ^/?2i  -ppr  Isi  =  nrjn 

Because  p  is  a  unit  vector,  we  recognize  that  p  •  p  =  p2  =  1 .  In  our  case,  we  can  express  the  vectors  in  the 
rotating  e-frame,  thus  the  quantity  (fix  p)  •  (fix  p)  may  be  computed  as  follows. 


er  e,  e„ 

-<j>c9s<j>  -\v  +  9^c(j>s9 

(Bxp)  = 

<i bsO  -(j>c9  v  +  9 

c(j>c9  c(j>s9  s</> 
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c(j>c9  (v  +  -  <j>s<j)s9 

(Qxp)-(Qxp)  =  if)1  c1 9s2  <j>  +  2<j)(\ >  +  c(j>c6s(j)s6  +  c2(j)s29[y  +  + 

c2</>c20(y  +  0)  - 2^(i>  +  9^c(j>c9s<j>s9  +  (j>2s29s2(f  +  <j>2c2<t> 

=  c2<j)(v  +  9 )  +(/)2 

This  result  is  of  course  the  same  as  that  achieved  using  the  specific  inertia  tensor  in  the  kinetic  energy 
formulation.  Using  the  equivalent  mass  f  tc  which  accounts  for  end-masses  and  the  tether  mass,  the  radius 
of  gyration  can  be  taken  as  that  of  a  very  thin  rod. 
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Appendix  D:  Operational  Limitations 


Choosing  a  proper  maximum  allowable  rms  control  current  is  important  to  ensure  feasible  EDT 
designs.  Permitting  too  much  current  can  cause  a  flexible  tether  to  curve  too  much,  tus  negating  the 
“dumbbell”  model  used  by  the  controller.  If  the  EDT  system  is  not  capable  of  driving  enough  average 
current  through  the  wire,  then  drag  will  overpower  the  electrodynamic  thrust  and  the  EDT  will  reenter  the 
atmosphere.  The  subject  of  this  appendix  is  to  determine  bounds  on  the  maximum  allowable  current  with 
which  we  may  constrain  the  optimal  control  problems  presented  in  Chapter  III  through  Chapter  V. 


Preventing  Reentry 


There  are  many  forces  acting  on  the  electrodynamic  tether  besides  gravity.  The  predominant  two 
are  the  drag  force  and,  when  the  electrodynamic  tether  is  actively  thrusting,  the  Lorenz  force  (see  Figure  1). 
If  the  cumulative  drag  force  is  greater  than  the  Lorenz  thrust,  then  the  EDT  will  reenter  the  atmosphere.  To 
determine  the  tether  system  design  requirements  to  compensate  for  atmospheric  drag  to  prevent  reentiy  we 

look  at  approximate  models  for  the  atmosphere  and  Earth’s  magnetic  field.  The  atmosphere  is  modeled  as 

* 

having  exponentially  decaying  air  density  with  a  scale  height  h  =  44  km  and  a  reference  altitude  and 
density  of  h0  =  250  km  and  po  =  6e  - 1 1  kg/nT  respectively.  In  a  circular  orbit  the  drag  force  will  always 


act  opposite  the  direction  of  motion  (i.e.  the  -et  direction).  The  drag  magnitude  on  the  two  endmasses  is 

(67) 


1 


-(*-*„/) 


1  (»' 


Fa ,  +Fal  =  -v-Poe  *•  (CdA  ♦  CdlA, )  =  ^  p{h)(CnA,  +  Cd2A2 ) 


2  v  r  j 


-(h-href) 

where  p{h)  =  poe  h  .  The  drag  force  on  the  tether  alone  is  approximated  using  a  constant  force 

distribution  along  the  tether.  This  approximation  is  adequate  for  this  calculation  when  the  COM  of  the 
system  is  close  to  the  midpoint  of  the  tether.  Justification  of  this  approximation  is  shown  as  follows. 
Letting  dA  =  d.ds  .  the  force  due  to  the  drag  on  the  tether  alone  is 
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'  tether 


L  2  ~S 

=  } dF  =  —Cdlv2 } pdA  =  —  Cdy p(h^dt  f  eh  ds  =  Cdtv~ p(h\dth  sinh 

2  2  -i/2 


2  h  ) 


Letting  z  =  — T  ,  an  expansion  of  the  hyperbolic  sine  function  is 
2h 

3 

z 

sinh(z)  =  z  +  —  +  h.o.t 
6 

If  we  are  willing  to  ignore  0(z3)  terms,  then  we  may  let  sinh(z)  =  z  and  write 


F,etHer=l-CdtV2p(h)dtL 

where  we  recognize  that  dtL  is  the  cross  sectional  area  of  the  tether.  Only  in  very  long  tethers 
( L  >  35  km  )  at  low  altitudes  ( h  <  250  km  )  could  the  0(z3)  terms  be  too  significant  to  ignore.  The  total 
aerodynamic  drag  force  on  the  entire  system  is  the  sum  of  the  drag  forces  on  both  end  masses  and  on  the 
tether  itself. 


D  =  |  v'p  (h)  (Fal  +  Fal  +  Fmher )  =  | 


p(h)(Cd]At+CdlA2  +  Cd,Ld,) 


(68) 


The  maximum  Lorenz  force  on  the  other  hand  depends  on  the  maximum  average  current  that  is 
permitted  through  the  wire  and,  for  the  nadir-pointing  tether  in  a  circular  orbit,  will  only  have  a  component 
of  thrust  in  the  positive  velocity  direction  et  perpendicular  to  the  normal  component  of  the  local  magnetic 
field. 

f  =  (7,AxB)'et 


For  a  circular  orbit  the  force  magnitude  is 


F  =  ImLB0 


(  r  \ 


V  r  J 


COS  I 


(69) 
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N 

where  the  circular  reference  orbit  is  taken  to  be  at  a  250  km  altitude  where  B  =  2.7 e  -  5  - .  The 

0  A-m 


magnetic  flux  density  component  normal  to  the  orbital  plane  decays  as  the  magnetic  inclination  i 
approaches  90  degrees,  and  thus  so  does  the  transverse  force  component  of  a  nadir-pointing  EDT.  For  a 
given  orbit  in  a  drag  environment,  the  average  Lorenz  force  must  be  greater  than  the  average  drag  or  the 
EDT  will  lose  altitude  quickly  and  burn  up  on  reentiy. 

It  is  desirable  to  graph  the  force  balance  altitudes  for  various  tether  designs  from  which  we  can 
perform  design  trade  studies.  Knowing  the  lowest  feasible  altitude  one  can  achieve  for  various 
electrodynamic  tether  designs,  one  can  determine  a  control  strategy  that  avoids  inadvertent  reentry.  The 
ratio  of  average  forces  given  in  Eq.  (68)  and  Eq.  (69)  is 


F 

D 


cos i 

i  u  ( r\  Hh~K) 

if\r)p°e  M+C^+W) 


Separating  the  reference  altitude  parameters  (ra,  Ba  and  p0 )  and  assuming  that  the  coefficients  of  drag  are 
all  approximately  2,  we  write  the  force  ratio  in  terms  of  the  radial  distance  or  equivalently  the  altitude. 


m  o 
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- 1 
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Vo  J  \ 

(70) 


where  A  =  Al  +  A2+  d,L  is  the  total  system  cross  sectional  area.  To  avoid  drag  induced  orbit  decay,  the 


average  available  Lorenz  force  must  overcome  the  average  drag  force,  i.e.  —  >  1 .  A  graph  depicting  the 

force  ratio  as  given  by  Eq.  (70)  is  shown  in  Figure  41  for  a  2  Amp  EDT  in  a  circular  equatorial  orbit  with 
the  force  balance  condition  depicted  as  a  vertical  dashed  line.  The  graph  clearly  shows  that  long  skinny 
tethers  permit  orbits  at  lower  altitudes  than  short  wide  tethers  which  are  more  susceptible  to  drag  and  do 
not  generate  as  much  thrust.  Varying  design  parameters  such  as  the  maximum  average  current  or  orbit 
parameters  such  as  the  inclination  will  shift  along  the  curves  according  to  Eq.  (70).  For  instance, 
increasing  the  orbit  inclination  reduces  average  normal  component  of  the  magnetic  flux  density.  This 
component  is  the  only  one  that  contributes  to  forward  thrust  in  a  nadir-pointing  EDT,  so  reducing  it 
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diminishes  the  system’s  boosting  capability  and  as  a  result  the  available  thrusting  force  to  overcome  drag  is 
reduced  thus  reducing  the  margin  to  the  force  balance  altitude  (or  forces  higher  force  balance  altitude). 
Increasing  the  maximum  allowable  control  current  increases  the  available  thrusting  force  and  shifts  the 
curves  right,  resulting  in  a  lower  force  balance  altitude  (or  increases  margin  to  force  balance  altitude).  This 
graph  is  use  till  for  design  trade  studies  to  determine  the  minimum  force  required  to  maintain  a  given  orbit. 
The  upper  bound  on  the  generated  force  is  limited  by  the  tether  quasi-equilibrium  curvature  described  in 
Appendix  D. 


Max  Lorenz  Force  to  Drag  Force  balance,  circ  orbit,  dt  =  4  mm 


Figure  41.  Maximum  Lorenz  Force  to  Drag  Force  Ratio  at  the  Magnetic  Equator 

It  is  should  be  recognized  that  there  are  other  considerations  that  factor  into  the  design 
requirements  such  as  ohmic  losses  in  the  conducting  wire,  cosine  losses  due  to  non-vertical  wire  orientation 
in  a  spinning  or  librating  tether,  and  even  non-operation  during  eclipse  times.  The  design  limits  presented 
here  represent  an  absolute  lower  bound  on  the  average  current,  therefore  the  peak  current  available  for  a 
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real  system  would  need  to  be  higher  than  this  average  value.  The  upper  bound  on  the  allowable  control 
current  is  driven  by  the  tether  length  and  diameter,  and  is  described  in  the  following  section. 


Validity  of  a  Straight  Tether  Model 


To  control  an  electrodynamic  tether,  it  is  important  to  understand  the  tether  shape.  Simpler 
control  laws  using  current  through  the  conducting  wire  are  available  when  we  assume  that  the  tether  lies 
along  a  straight  line  between  the  two  end-masses.  To  justify  this  approximation,  we  need  an  adequate 
shape  model  with  which  we  can  determine  the  tether  constraints  that  maintain  a  relatively  straight  line.  For 
a  given  orbit,  the  EDT’s  maximum  current,  length  and  diameter  all  factor  into  the  tether  shape  and 
vibration  dynamics.  These  parameters  must  be  considered  to  ensure  feasible  control  solutions.  An 
approximate  model  using  spectral  separation  developed  by  Von  Flotow  (Ref  54  pp  76-77)  will  serve  this 
purpose.  An  outline  of  these  approximations  is  given  here. 

Because  the  dynamics  of  the  flexible  tether  experiences  fast  motion  (longitudinal  vibration  along 
the  tether)  and  slow  motion  (lateral  vibration),  we  may  view  the  slow  dynamics  as  forming  an  equilibrium 
with  respect  to  the  fast  dynamics.  A  quasi-equilibrium  state  may  be  reached  in  slow  time  when  the  lateral 
force  distribution  along  a  vertical  wire  is  balanced  by  gravity  gradient  and  tension  forces.  Viewing  the 
quasi-equilibrium  in  this  manner  permits  us  to  determine  an  instantaneous  tether  shape  using  the  following 
equation. 


=  0 


where  S  is  the  distance  along  the  strained  arc-length  of  the  tether  subject  to  tension,  r  ,  and  external  forces 
per  unit  length,  fext ,  such  as  drag  and  Lorenz  forces  (Figure  42). 
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Figure  42.  Tether  Curvature  Due  To  Lorenz  and  Aerodynamic  Force  Distribution 

The  vector  position  of  the  tether  is  R  ,  thus  the  equilibrium  radius  of  curvature  is  its  magnitude  R  .  The 
tether  does  not  curve  very  much  and  we  can  approximate  it  as  that  of  a  circular  arc  of  radius  R  .  If  we 
assume  that  the  tether  mass  is  much  less  than  the  end-body  masses,  then  we  can  assume  the  tension  is 
independent  of  S  .  At  equilibrium  conditions  with  these  assumptions,  the  tether  radius  of  curvature  is 

T  . 

R  =  —  where  ft  is  the  total  external  force  density  component  in  the  et  direction  perpendicular  to  the  line 

ft 

between  to,  and  m2 .  This  lateral  force  density  is  measured  in  units  of  force  per  unit  length.  Tether 
curvature  in  the  er  -  et  plane  is  depicted  in  Figure  43  with  (j)  =  0 . 

The  curve  angle  <//  may  be  written 

2¥  =  ^V'  =  ^-  (71) 

R  2  t 
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Figure  43.  Curved  Tether  Geometry 


The  tether  shape  is  determined  by  its  tension,  length  and  transverse  component  of  external  forces  acting  on 
it.  To  justify  a  straight-line  tether  assumption  in  the  dynamic  equations  and  subsequent  control  laws,  we 
cannot  allow  the  tether  to  curve  appreciably  under  the  distributed  loads  along  the  tether  to  the  extent  that 
the  resultant  force  is  significantly  smaller  than  that  of  a  straight  tether.  From  Figure  43  we  see  that  the  arc 
portion  Ry/  is  slightly  longer  than  the  vertical  component  of  the  wire  that  effectively  produces  thrust  in  the 

p±  direction,  R  sin  i//  .  This  difference  results  in  the  straight-line  model  error  given  by 


eslm=R(siny/-y/)  =  R  y/-—-y/  +  h.o.t. 

6 


Ri//3 

~ 


where  we  have  made  use  of  a  Taylor  expansion  for  y/  about  zero  and  ignored  terms  higher  than  order 
three.  If  we  desire  the  straight  tether  thrust  approximation  to  remain  within  95%  that  of  the  actual  curved 
tether,  we  require  that 

R  sin  y/>  0.95*  Ry/  =>  sin  y/  >  Q.95y/ 


therefore, 
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swlim^lim--^  =  0.95^lim 

6 

*'i-6(1-95)-0-3 

V'lim  =  ±0.548 


To  obtain  the  force  per  unit  length  /'  we  include  electrodynamic  and  atmospheric  effects,  fte  and  fa .  The 
maximum  electrodynamic  force  per  unit  length  for  a  given  7max  is 


fe  =■ 


7  LxB 


7 max  Bn 


where  /?.  is  the  en  component  in  the  rotating  frame  of  the  Earth’s  magnetic  flux  density.  At  250  km,  this 
magnetic  flux  density  is  about  2.7  e-5  N/(A-m)  (rounded  average  from  International  Geomagnetic 
Reference  Field). 

The  maximum  aerodynamic  force  per  unit  length  of  the  wire  in  nadir-pointing  equilibrium  and  circular 
orbit  is 


4  =  v2p(h)cosad,  0  =  ^~p(h)d, 


so  the  total  lateral  force  density  is 

f,=  f,e+ f, a  =  “4a x4  +  ^p{h)d, 

At  250  km  where  the  atmospheric  density  is  about  6e-l  1  kg/m3,  this  force  density  is  about  1.8e-5  N/m. 
Using  t/xlimto  determine  f]  using  Eq.  (71)  then  substituting  into  Eq.  (72)  yields  a  maximum  allowable 
current  that  permits  the  straight  tether  line  assumption. 


(72) 


1 


4iax  . , 

Bn  V 


2-^-^p(h  k 

4  Cl 


For  a  2  km  tether  in  a  250  km  orbit  with  a  tension  at  nadir-pointing  equilibrium  of  about  .7  N,  this 
maximum  current  equates  to  about  13  Amps.  A  tether  carrying  13  Amps  of  current  will  have  a  radius  of 
curvature  and  mid-point  sag  distance  given  respectively  by 
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.7 

13-(2.7e-5) 


=  1920  m 


8 


—  =  260  m 
8  R 


Von  Flotow’s  paper  provides  a  formula  for  the  period  of  the  first  lateral  mode  of  vibration  as 


Jfl  / 

where  8t  is  the  mass  density  of  the  tether.  For  a  uniform  density  tether,  this  is  simply  8t=  %  and  the 
period  of  lateral  vibration  is 

For  a  400  kg  tether  system  described  above  (10%  of  which  is  tether  mass),  this  period  would  equate  to  676 
seconds.  Flaving  bounded  the  maximum  allowable  control  current,  we  may  now  pose  optimal  control 
problems  that  can  achieve  feasible  solutions. 
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Appendix  E:  Tether  Tension  Curves 


There  are  two  sets  of  graphs  in  this  Appendix.  The  first  set  represents  solutions  to  equation  (40) 
for  various  altitudes,  mass  ratios,  ballistic  coefficient  ratios  and  tether  lengths.  In  all  cases,  the  numerical 
solutions  use  a  400  kg  system  with  a  tether  thickness  of  5  mm.  The  second  set  of  graphs  show  the  in-plane 
equilibrium  conditions  at  which  the  tether  tension  is  zero.  Equilibrium  conditions  above  the  curves  are 
where  the  tether  is  in  tension,  while  regions  under  the  curves  represent  non-tension  conditions  where  the 
tether  would  be  slack.  The  slack  conditions  must  be  avoided  when  determining  the  EDT  control  strategy  in 
which  case  we  may  bound  the  Estate  to  meet  the  constraints  established  by  Eq.  (39)  for  all  libration  angles 
and  not  just  at  equilibrium. 


114 


Tension  [N]  Tension  [N] 


Tension  of  tethers  in  equilibrium  at  150  km  circ  orbit  for2/3<  (B1*/B2*)  <1 


Equilibrium  Theta  [deg] 


Tension  of  tethers  in  equilibrium  at  150  km  circ  orbit  for  B1*/B2*  =  3/2 
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Tension  [N]  Tension  [N] 


Tension  of  tethers  in  equilibrium  at  150  km  circ  orbit  for2/3<  (B17B2*)  <1 


Equilibrium  Theta  [deg] 


Tension  of  tethers  in  equilibrium  at  150  km  circ  orbit  for  B17B2*  =  3/2 


Equilibrium  Theta  [deg] 
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Tension  [N]  Tension  [N] 


Tension  of  tethers  in  equilibrium  at  200  km  circ  orbit  for2/3<  (B1*/B2*)  <3/2 


Equilibrium  Theta  [deg] 


Tension  of  tethers  in  equilibrium  at  200  km  circ  orbit  for2/3<  (B1*/B2*)  <3/2 


Equilibrium  Theta  [deg] 
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Tension  [N]  Tension  [N] 


Tension  of  tethers  in  equilibrium  at  250  km  circ  orbit  for2/3<  (B1*/B2*)  <3/2 


Equilibrium  Theta  [deg] 


Tension  of  tethers  in  equilibrium  at  250  km  circ  orbit  for2/3<  (B1*/B2*)  <3/2 


Equilibrium  Theta  [deg] 
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Tension  [N]  Tension  [N] 


Tension  of  tethers  in  equilibrium  at  300  km  circ  orbit  for2/3<  (B1*/B2*)  <3/2 


Tension  of  tethers  in  equilibrium  at  300  km  circ  orbit  for2/3<  (B1*/B2*)  <3/2 


Equilibrium  Theta  [deg] 
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Tether  Length  [km]  Tether  Length  [km] 


In-plane  Lib  Theta  @  zero  tension,  150  km  alt  circ  orbit 


In-plane  Lib  Theta  @  zero  tension,  150  km  alt  circ  orbit 
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Tether  Length  [km]  Tether  Length  [km] 


In-plane  Lib  Theta  @  zero  tension,  150  km  alt  circ  orbit 


Zero  Tension  Theta  [deg] 


In-plane  Lib  Theta  @  zero  tension,  250  km  alt  circ  orbit 


Zero  Tension  Theta  [deg] 
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Tether  Length  [km]  Tether  Length  [km] 


10 


In-plane  Lib  Theta  @  zero  tension,  250  km  alt  circ  orbit,  B17B2*  =  2/3 
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Appendix  F:  Scaling 


Because  of  the  nature  of  the  optimal  control  problems  posed  in  this  paper,  there  are  a  number  of 
states  that  can  be  very  large  or  very  small  and  ones  that  change  over  small  and  large  time  durations. 
Furthermore,  it  is  important  to  scale  the  problem  parameters,  not  only  to  assist  in  derivations  of  equations 
of  motion,  but  also  to  condition  a  problem  to  achieve  an  accurate  numerical  solution.  This  appendix 
addresses  the  scaling  that  is  required  to  achieve  both  goals.  First,  scaling  the  time  variable  to  derive 
averaged  state  equations  of  motion,  and  second  to  scale  all  the  problem  parameters  for  well  conditioned 
numerical  solutions. 


Scaling  the  Time  Variable  for  Derivation  of  the  Averaged  State 
Equations  of  Motion 

The  control  problems  posed  in  this  paper  include  states  that  undergo  small  rapid  changes  over 
short  time  scales,  but  on  average  change  more  significantly  over  long  time  scales.  For  maneuvers 
spanning  long  time  scales,  it  is  advantageous  to  average  out  the  small  fluctuations  occurring  over  the  short 
time  scales  and  only  consider  the  long  term  behavior  of  the  average  states.  This  is  achieved  through  the 
method  of  averaging  offered  by  perturbation  theory.  Averaging  a  state  over  a  2 n  period  requires 
integration  of  the  instantaneous  state  with  respect  to  a  time  variable.  It  is  important  to  recognize  which 
terms  change  rapidly  and  must  be  integrated  and  which  terms  change  so  slowly  that  they  may  be  considered 
constant  over  a  single  period.  To  assist  in  this  clarification,  two  time  scales  are  employed  to  identify 
parameters  that  change  slowly  and  ones  that  change  quickly. 

The  true  anomaly,  v  ,  is  related  to  the  clock  time  t  through  the  equation 

v  =  nt 

where  n  is  the  mean  motion  of  the  orbit.  This  variable  v  is  be  referred  to  as  a  “fast”  time 
variable  which  changes  rapidly  on  a  “short”  time  scale,  e.g.  over  a  single  orbit.  On  the  other  hand,  the 
variable  T  is  related  to  the  clock  time  t  by 
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T  =  £t  =  £V/ 

where  £  is  a  small  parameter  such  that  s  <K  1 .  In  this  paper,  a  useful  scaling  factor  is  the  non-dimensional 
electrodynamic  torque  defined  in  Eq.  (43).  The  variable  T  captures  the  dynamics  which  vary  slowly  and  is 
referred  to  as  the  “slow”  time  variable.  Because  T  changes  slowly  with  time,  it  represents  the  dynamics 
which  take  place  on  the  “long”  time  scale,  i.e.  over  many  orbits.  These  variables  may  be  treated  as 
independent  variables  so  long  as  s  <K  1 .  Although  two  time  scales  are  employed  in  this  paper,  multiple 
time  scales  may  also  be  used  in  a  similar  fashion  (e.g.  including  a  rotating  tilted  magnetic  dipole  moment 
would  introduce  a  medium  time  scale  variable). 


Scaling  the  OCP  for  Well-Conditioned  Numerical  Solutions 

When  solving  an  optimal  control  problem  (OCP)  using  numerical  methods,  it  is  critical  to 
condition  the  input  parameters  to  achieve  accurate  results  with  faster  CPU  run  times.  Scaling  is  essential  to 
writing  a  well  conditioned  problem  for  the  computer.  This  section  draws  from  discussions  given  in  Ref  41, 
pg  31-36,  and  Ref  47  the  salient  points  being  repeated  here  for  convenience.  Recall  the  unsealed  state  and 

control  vectors,  \{T)  =  [a,h,k,i,z]T  and  u(7’)  =  [nl,u1,u3,u4,u5]T  respectively.  Numerically,  it  is 
advantageous  to  scale  each  element  in  the  state  and  control  vectors  for  computational  efficiency,  i.e. 

=  [a,h,k,i  ,zj 

[  R  5  ^ 2  5  ^3  5  ^4  dd  ] 


where  A,H,K,I ,Z,U1,U2,U3,U4,  and  U5  are  arbitrary  designer  units.  For  the  problems  posed  in  this 
paper,  designer  units  chosen  to  make  the  scaled  states  and  controls  roughly  order  one  worked  well.  Time  is 
scaled  in  a  similar  fashion  expressed  using  a  designer  clock  time  unit,  ts . 

t 


u(j)  = 


u ,  u ,  ux  UA 


Ur 


u2  "3  UA  **5 

c/1’c/2,c/3,c/4,c/5 


7\T)= 


a  h  k  i  z 

H’H’-R’T’I 
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For  the  work  presented  here,  the  designer  clock  time  unit  was  chosen  such  that  tf  is  order  one.  As  an 
example,  the  OCP  presented  in  Eq.  (17)  is  rewritten  in  scaled  form  to  input  into  the  numerical  optimizer. 
Minimize  Cost:  J  =  t,  (73) 

Subject  to: 


t‘  ts  ts  ts  lf(x(7’),u(r)) 


a’h’k’i’z 


i(r)  =  diag 
e0  (X  (^0  ))  =  (0’  zo] 

&(u(r))  =  £,-/’<<) 

where  /  is  the  maximum  allowable  rms  control  current.  The  scaled  box  constraints  are  enforced  as  well. 


x,  <x(t  )<x„ 
u,  <u(T)<u, 

_  _  _  (74) 


Note  that  the  OCP  presented  here  is  mathematically  identical  to  the  unsealed  one  in  Eq.  (17). 

Also,  recognize  that  the  dynamic  constraint  vector  f  has  elements  containing  unsealed  states.  The  output, 
however,  is  scaled  to  be  compatible  with  the  other  scaled  constraints.  Using  the  fact  that  a  generic  unsealed 
state  and  time  are  related  to  their  scaled  counterparts  by  x  =  xX  and  1  =  tts ,  each  individual  time  derivative 
is  scaled  as 


dx  dx  dx  dt  .  ts 
dt  dt  dx  dT  X 


The  scaled  OCP  in  the  form  of  Eqs.  (73)  and  (74)  were  used  as  the  input  for  the  optimizer  (DIDO)  that 
yielded  the  solutions  in  this  paper. 
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Appendix  G:  Derivation  of  Averaged  Orbital  Element 
Equations  of  Motion 


This  appendix  provides  a  derivation  of  secular  equations  of  motion  using  a  mixed  set  of  classical  and 
equinoctial  coordinates  and  the  method  of  averaging.  To  determine  the  secular  change  in  a  given  state  xt , 
we  start  with  the  perturbation  equations  of  motion  given  in  Eq.(9)  and  use  the  approximation 

dt  »  — (l-2ecos v)dv  then  integrate  over  N periods  as  follows. 
n 

rt.  1  f2/rV  dx .  /  , 

Ax ,  =  \f  dx.  — L(l-2ecosvW 

•''»  n  ■’°  dt  v 


Because  the  orbits  considered  are  nearly  circular,  eccentricity  is  very  small  and  the  argument  of  perigee 
is  ill  defined.  Therefore  two  equinoctial  coordinates  defined  as  h  =  e  s  i  n  o>  and  k  =  <?cos  co  are  better  suited 
for  this  orbit  type.  Thus,  changes  in  semi-major  axis,  inclination  and  right  ascension  of  the  ascending  node 
are  approximated  as 


where  C  =  - 


B*fip(a )  Cl 

-  represents  the  thrust  per  unit  current  and  the  drag  rate  is  D  = - -  .  Note  that  —  is 


dimensionless. 


The  only  control  that  will  yield  non-zero  solutions  after  integrating  the  above  equations  is  a  periodic 
current.  The  control  current  may  be  expressed  as  the  sum  of  the  periodic  functions  that  produce  secular 
changes  to  the  states,  therefore  we  use  the  first  five  terms  of  a  Fourier  series  shown  in  Eq.  (1).  After 
integration  we  obtain  the  secular  changes  to  three  of  the  five  states  that  change  on  a  long  time  scale. 
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A  a  «  [2C/macosj'(M,  +  ;/2e)-2Z)] 


2nN 


A  l  a 

AQ 


— C7_  sin  i 


ux  k2-h2 


hk 


4e2  ’ 4  2e2  5 


-  -C/„ 


hk  k2-h 2 

2^U4 


2e 


4e- 


2nN 


2jtN_ 

y  « 


(75) 


The  time  derivatives  of  the  equinoctial  coordinates  may  be  calculated  as  follows. 


h  =  ecosoro  +  esin® 

«  Cl  cos  ;  [cos  [v  +  ®)(/?cos  v  +  £sinv)(l  +  2ecosv)cos®  +  (2  -ecos  v)(l  +  3ecosv')sinv'cos® 

+  (2cosv +  esin2  v)(l  +  3ecosv')sinol-^--  l  +  fl  +  ^lecosK  (2-ecosv')cos®  (76) 

v  J  a  [  2e  [  ^  i  J 


+ 

( 

COS  V  +  6?  + 1 

ecos2  v 

sin 

V 

l  h  j 

J 

!  i 

Carrying  out  the  multiplications,  eliminating  second  and  higher  order  eccentricity  terms,  then 
substituting  in  h  and  k  ,  we  write 


h  «  Cl  cos; 
2h 


k  h  T  (  5 

— cost” - sini/  +  2&cos2  v'-/7sin2v  (/?cosv'  +  £sinv')  +  2sini/ +  — esin2v 

e  e  P  M  2 


+  —  cos  v  +  h  +  5 h  cos  v  >  - 
e  a 


2D 


-  +  — -sin  2k-  1 H — —  —  cos'  v  sin  v 
e  2 h  \  h  )2  ]  e 


+ 

( 

COS  V  +  6?  +  | 

fi+4] 

ecos2  v 

V 

l  h  ) 

1  J 

(77) 


k  h 

recognizing  that  cos(v  +  co)  =  cosKeosfu-sin vsin&>  =  —  cosk - sinv  and  e2=h2+k2. 

e  e 

Integrating  with  respect  to  the  true  anomaly  from  0  to  2nN  ,  we  find  the  change  in  the  average  h  state. 


Ah  «  <  Cl  cos  i 


3  h\  (h)  (k 
u,  ]  —  |  + ]  —  |  +  u,  |  —  [  +  u, 
e 


1  2  r2 


rh  M2^ 
v 4  +  2e:  j 


+  u , 


k  ( k2-h2)k 

4+  4e2 

v  j 


-q1  +  4 V4— (78) 
a  (  h  )  \  n 


We  obtain  the  k  state  dynamics  in  a  similar  manner. 
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k  =  -e  sin  + ecos® 

«  Cl cos;[-cos(v  +  ®)(/?cosv  +  £sinv)(l  +  2ecosv)sin®-(2-ecosv)(l  +  3ecosv)sinvsin® 


(2cosv  +  esin2 v)(l  +  3ecosv)cos®j--^J  SmV  1  +  j  1  +  7  jecosv  (2-ecosv)sin®  (79) 

Cl  I  2^  ,  (  h  )  y 


1  u  1 

+  cosv  +  e  +  1  +  —  ecos  v  cos co 

,  V  h  . 


Rewriting  Eq.  (79)  using  definitions  of  h  and  k  ,  we  obtain 


k  «  Cl  cos;  j  cosv-  — sinv  +  2£cos2  v-/7sin2vj(/?cosv  +  £sinv)-^2sinv  +  -^-esin2vj  — 

2k  ,  ,,  ,  ]  2 D\  ( shiv  a  .  „  ( ,  a\e  ,  •  V* 

+  —  cosv  +k  +  5jfccos~  v  1 - <-  - H - -sin2v-  1  +  — r  —  cos’vsmv  — 

e  \  a  \  y  e  2h  V  h  )  2  )e 

(  f,  a)  2  Ul 

+  cosv  +  e+  1  +  —  ecos  v  —> 

\h  e 


Integrating  with  respect  to  the  true  anomaly  from  0  to  2nN  ,  we  find  the  change  in  the  k  state. 


A,  I  ■  r (k\  f-h)  fk  h2k)  f  h  [h2-k2)h\  Df  a\  2nN 

Ak « j cim cos ;  "  - 7  1+*Tnr 
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Appendix  H:  Propagation  of  Libration 


When  propagating  a  variable  such  as  libration  that  exhibits  rapidly  changing  behavior  on  short 
time  scales  but  also  exhibits  slowly  changing  behavior  over  longer  time  scales,  it  is  necessary  to  use  stiff 
ordinary  differential  equation  (ode)  solvers.  These  numerical  solvers  are  subject  to  errors  which  can  grow 
as  solutions  are  propagated  over  a  long  interval.  To  test  the  ode  solver  in  a  problem  relevant  to  this 
research,  the  following  homogeneous  ordinary  differential  equation  was  propagated. 

0(v)  +  40(v)  =  O 

(oJ 

Mi =o 


revs 


Figure  44.  Matlab  ode23  Solution  to  Homogeneous  Equation  vs.  Exact  Solution 

The  propagated  solution  to  Eq.  (81)  is  plotted  in  Figure  44  along  with  the  exact  solution  to  this  equation 
using  Matlab’s  ode23t  (stiff  solver).  Notice  that  over  the  course  of  time  the  propagated  solution’s  phase 
slowly  drifts  from  the  exact  solution  (i.e.  (/)  =  (j)m  cos2v  )  where  1  rev  =  2 n  rad.  The  shift  is  not  due  to  real 
perturbations  since  this  is  an  exact  homogeneous  solution,  but  rather  due  to  numerical  round  off  errors 
which  must  be  addressed  when  propagating  control  solutions.  When  using  a  propagator  to  model  a 
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nonlinear  plant  response  to  a  control  input,  the  controller  could  become  out  of  phase  with  the  propagated 
states  due  solely  to  numerical  errors.  To  compensate  for  this  small  artificial  phase  shift  when  propagating 
control  solutions,  we  provide  the  controller  given  by  Eq.  (1)  with  a  slightly  phase  shifted  true  anomaly 
input,  v c ,  at  each  instant.  This  correction  has  the  effect  of  slightly  raising  or  lowering  the  frequency  of  the 
periodic  controller  to  match  the  frequency  of  the  numerically  propagated  homogeneous  solution.  The 
modified  true  anomaly  is  designed  to  be  in  phase  with  the  propagator  to  simulate  a  real  plant  that  is 
unaffected  by  round  off  error  and  is  defined  by 

f  Av 

vc=v  1--^  (82) 

l  a  y) 

A  v 

where  - —  is  the  change  in  the  propagated  independent  variable  v  with  respect  to  the  change  in  the  true 

A  v 

independent  variable  v  .  The  easiest  way  to  obtain  A  v  is  to  determine  the  difference  between  the  v  (in 
revs)  corresponding  to  the  final  peak  of  the  propagated  homogeneous  solution  after  A  v  =  N  revs, 

( v„ )  and  the  peak  of  the  exact  solution  (always  corresponding  to  either  whole  or  half  revs).  The 

difference  A  v p  is  positive  when  the  propagated  peak  lies  to  the  right  of  the  exact  solution  peak  (i.e.  lags  the 
exact  solution)  and  is  negative  when  it  lies  to  the  left  (i.e.  leads  exact  solution),  written  as 

Ay  (vp)peak~N 
Av  N 

The  value  is  non-dimensional  and  will  shift  the  control  input  variable  vc  throughout  the  trajectory  to 
maintain  phase  with  the  propagated  phase  vp  (i.e.  control  frequency  is  matched  with  the  “natural” 

frequency  as  determined  by  the  numerical  propagator).  Note  that  at  the  final  time  after  N  revs  the 
controller  is  in  phase  with  the  propagated  trajectory,  i.e. 

(  Av 

v=N  1 - p-  =N  -  Av 

N  p 

v  JV  J 

The  controller  is  then  rewritten  as 

/  =h,  (T)  +  u2  (J)eosvc  +u3  (J)sinvc  +u4  (J)eos2vc  +u5  (r)sin2v,  (83) 
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It  must  be  emphasized  that  this  controller  is  only  used  for  propagation  to  compensate  for  numerical  errors 
to  achieve  more  accurate  comparisons.  This  scheme  is  not  necessary  when  applied  to  a  real  world  design, 
although  some  variation  of  this  method  would  be  useful  for  real  perturbations  originating  from  other 
sources. 
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Appendix  I:  Reference  Synopses 


Reference  Author 

JGCD  Vol  26,  #5  Sep-Oct  2003  Tragesser  &  San  Orbital  Maneuvering  with  Electrodynamic  Tethers 

Titlted  Dipole  Magnetic  field.  Perturbation  equs  (a,e,I,...)  w/  expansion  in  e  appx  to  get  secular 
changes.  Change  shows  how  various  current  laws  change  a,  e,  I,  etc.  Ex  I(nu)  =  cos(nu)  changes  e. 
EOM  assumes  tether  aligned  w/  vertical.  Rigid  rod  model.  Changes  indep  variable  from  t  to  nu. 
Maneuvers  are  NOT  optimal. 

JGCD  Vol  28  #2  P  Williams  Optimal  Orbital  Transfer  with  Electrodynamic  Tether 

Change  orbit  by  modulating  current  in  wire.  Takes  into  account  librations,  both  in  plane  and  out  of 
plane.  Compares  w/  similar  maneuvers  using  a  hypothetical  non-librating  tether.  Has  useful  reduced 
mass.  No  atmosphere  considered.  Example  problem  uses  a  500  km  orbit. 


Title  and  Synopsis 


AIAA-2001-1 139 


Life  Extension  and  Orbit  Maneuvering  Strategies  for  Small  Satellites  in  LEO  Using  Electrodynamic 
Tethers. 

Has  good  sample  parameters  w/  which  to  frame  the  problem.  Addresses  micometeor  problem  and 
pourous  tape.  Shows  graphs  of  responses  to  various  tether  configurations  and  conditions.  No  EOMs 


AIAA- 1999-284 1-933  Gilchrist  Space  Electrodynamic  Tether  Propulsion  Technology:  System  Considerations  and  Future  Plans 

Good  plots  of  bare  wire  EDTs.  ProSEDS  and  TSS-1R  missions. 

AIAA-2000-440-65 1  Estes  Performance  and  Dynamics  of  an  Electrodynamic  Tether 

Discusses  advantages  of  a  bare  EDT  in  collecting  electrons.  Reviews  boosting  and  deboosting 
applications.  Gives  system  performance  variation  vs.  key  parameters.  Discusses  performance  of  bare 
EDTs  for  boost  or  de-boost  applications 

AIAA-200 1-3980- 120  Van  Noord  Electrodynamic  Tether  Optimization  for  the  STEP  AIRSEDS  mission 

Mostly  about  design,  survivability  and  manufacture  of  EDT.  No  EOM 

ataa  onm  \Aicsn  o  *  i  Evolution,  Technology  and  Direction  for  the  ASTOWSTEP  AIRSEDS  Electrodynamic  Tether 

AIAA-2003-143-567  Santangelo 

Mission 

Brief  format,  not  paper 

AIAA-2004-3501  -989  Vaughn  Review  of  the  PROSEDS  Electrodynamic  Tether  Mission  Development 

Mission  never  got  off  the  ground,  but  good  passdown  for  future  EDT  mission  planners. 

AIAA-2004-5309-275  Palaez  Self-Balanced  Electrodynamic  Tethers 

Inclined  orbits  produce  instabilities  on  EDTs.  Inert  tethers  are  fine  (for  circular  case).  Zeroizing  the 
Lorentz  torque  is  addressed  which  eliminates  instability.  Mass  distribution  is  critical,  else  damping  or 
control  is  needed.  Attitude  dynamics  and  Mag  field  model.  Rigid  rod  assumed. 


AIAA-2004-531 3-1 57  Watanabe  An  Application  of  Input  Shaping  for  Electrodynamic  Tether  Systems 

Input  shaping  to  reduce  vibrations  and  librations  on  an  EDT  being  propelled.  Mag  field  model 
includes  massive  flexible  lumped  mass  tether.  Considers  one  flexible  mode  of  vibration.  Uses  Kanes 
equation  (refs  4&7)  although  not  supplied.  Discusses  how  to  do  bang  bang  EDT  thrust  control  by 
stepping  in  a  controlled  way  that  reduces  vibes. 
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Reference 

Author 

Title  and  Synopsis 

Compares  3  different  electron  emission  methods  for  different  missions 

AIAA-2005-4545-358 

Bonometti 

Free  Reboost  Electrodynamic  Tether  on  the  International  Space  Station 

Proposes  tether  flywheel  design  to  reboost  ISS.  Has  some  charts  w/  ISS  data.  Altitude  history  1998- 
2004.  No  eom 

AIAA  2002-4641 

Tragesser 

Orbital  Design  of  Earth-Oriented  Tethered  Satellite  Formations 

Looks  at  dynamics  of  multiple  tethered  satellites  (3-D).  Investigates  stability.  Flexible  lumped  mass 
tether.  Starts  w/  rigid  model,  then  moves  to  flexible  elastic  one  uses  eoms  and  stability  analysis  from 
Hughes  book.  Example  formation  uses  circular  orbit. 

AIAA  5479-983 

Kumar 

Review  on  Dynamics  and  Control  of  Non  EDT  Satellite  Systems 

Light  on  equations  but  very  thorough  presentation  of  varius  work  in  dynamics  and  control  being  done 
by  different  researchers.  275  refs! 

AIAA  6473-285 

Lorenzini 

Libration  control  of  EDTs  in  Inclined  Orbit 

EOMS  derived  from  Lagrangian 

AIAA  6685-973 

Pelaez 

Dynamic  Stability  of  EDTs  in  Inclined  Elliptical  Orbits 

Elliptic  orbits  yield  periodic  solutions  not  equilibrium  positions.  EDT  control  can  be  employed  to 
manage  the  instability  for  small  eccentricity  orbits  (e<.35).  3  stages.  1.  Analyze  stability  of  elliptical 
orbit  inert  tether.  2.  Consider  electrodynamic  forces.  3.  Compare  w/  circular  case. 

AIAA  5077-785 

Pelaez 

Periodic  Solutions  in  EDTs  on  Inclined  Orbits 

Periodic  solutions  obtained  using  eigval  of  monodromy  matrices  compared  with  propagations  based 
on  Poincare  method  in  both  ep  (mag  to  grav  torque  ratio)  and  then  I  (inclination).  Even  compares  w/ 
linearized  solution.  Model  is  a  rigid  rod,  dominant  end  mass,  constant  tehter  current.  Periodic 
solution  exhibits  frequency  entrainment  phenomenon  (periodic  sol'n  has  same  period  or  integer 
multiple  as  forcing  terms) 

AIAA  17499-711 

Williams 

Libration  Control  of  Tethered  Satellites  in  Elliptical  Orbits 

Non-EDTs.  Control  via  tether  tension  (length  variation).  Stability  analyzed  via  floquet  theory. 

AIAA  4992-661 

Palaez 

Two  Bar  Model  for  Dynamics  and  Stability  of  EDTs 

Looks  at  2  rigid  bar  model  for  tether  to  analyze  the  impact  of  lateral  dynamics  on  stability.  Assumes 
Massive  s/c,  circ  orbit,  inclined,  but  only  constant  current,  no  control.  2  cases,  1.  Continuous 
conductive  wire,  2.  Part  discontinuous. 

AIAA  1990-1197 

Vadali 

Feedback  Control  of  Tethered  Satellites  Using  Lyapunov  Stability  Theory 
(also  hardcopy  of  Journal  0731-5090  Vol  14  #4  (729-735)) 

Has  Dynamic  Equations  useful  for  stationkeeping,  but  mostly  concentrates  on  deployment  and 
retrieval  of  tethers.  No  atm  drag.  Lyapunav  function  provided  for  in  plane  theta  and  L  eom. 
Coordinate  xform  is  presented  which  partially  decouples  theta  and  phi  dynamics. 

AIAA  9166-681 

Kojima 

Non-linear  Control  of  Librational  Motion  of  Tethered  Satellites  in  Elliptic  Orbits 

Controls  w/  thrusters  at  endmasses.  Assumes  no  aerodrag,  no  elasticity  of  tether,  only  in-plane 
libration  considered.  Mother  w/  2  subsatellites  connected  w/  massless  rigid  rod  tethers. 

AIAA  2002  4045 

Hoyt 

Stabilization  of  EDTs 

Uses  sensors  to  provide  feedback  control  varying  the  current  to  stabilize  the  EDT.  No  eoms,  but 
show  output  plots. 

AIAA  2000-322 

Lorenzini 

An  Overview  of  EDTs 

Good  overview  of  actual  deployed  systems  (TSS,  ProSEDS,  etc)  their  history,  what  we've  learned, 
what  we  can  do  in  future  missions 

AIAA  1990-656 

Matteis 

Dynamics  of  a  Tethered  Satellite  Subjected  to  Aerodynamic  Forces 
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AIAA  1990-1198 


AIAA  21768-992 


AIAA  1991-532 


AIAA  10546-681 


AIAA  11822-945 


AIAA  1990-1195 


AIAA  1991-1002 


AIAA-6934-481 
JGCD  2005  vol28#3  541 


AIAA  13956-480 


AIAA  2003-5781 


AIAA  4057-325 


Includes  aero  forces  on  the  subsatellite.  Aerodynamic  forces  play  a  role  in  determining  stability  of 
system  equilibrium.  Circ  orbit,  equatorial  plane,  rotating  atmosphere.  2  approaches  1.  Linearize  the 
eom.  2.  Propagate  non-linear  equations.  EOMS  assume  a  dominant  mass  (shuttle),  circ  orbit,  flexible 
tether.  Contains  excellent  aero  tether  refs. 


Von  Flotow  Insights  and  approximations  in  Dynamic  Analysis  of  Spacecraft  Tethers 

Discusses  vibrational  motion.  Equilibrium  shape  of  tether  is  slightly  sagged  from  straight  line. 
Includes  stretch  and  flexibility  in  tether.  Weak  instabilities.  Concludes  passive  damping  has 
inconclusive  effects. 

Yu  Periodic  Motion  in  the  Tethered  Satellite  System 

Motion  and  control  of  Mother/daughter  pair  in  circular  orbit,  then  elliptical  orbit.  Control  via  length 
rate  of  tether.  EOM  include  length  rate  and  tension  but  in  plane  motion  only.  Below  critical 
eccentricity  (e<.3)  motion  is  stable.  Limit  cycle. 

Matteis  Dynamics  of  a  Tethered  Satellite  in  Elliptical,  Non-Equatorial  Orbits 

(also  a  JGCD  article  0731-5090  1992  vol  15#3  (621-626).  Depicts  eoms  with  states  expressed  in 
Cartesian  coords.  2  cases.  1.  Equatorial  eccentric  orbits.  2.  Inclined  orbits.  Includes  aero  forces  and 
mentions  peak  natural  freqs  and  driving  aero  force  freq. 


Sanyal  Stability  and  Stabilization  of  Relative  Equilibria  of  Dumbell  Bodies  in  Central  Gravity 

Two  identical  masses,  rigid  rod,  massless  link.  5-DOF  eom.  Uses  Lagrangian  eom  approach  and 
then  Routh  reduction  to  eliminate  a  DOF. 

Somenzi  Linear  Stability  Analysis  of  EDTs 

Assume  circ  orbit,  inextensible  tether  2  pt  endmass.  Electrodynamic  forces  cause  coupling  of  cable 
oscillation.  Constant  current.  Bending  tether  under  current  load  included.  Lat  and  long  modes  of 
vibration.  Separates  lateral  modes  of  vib  from  librations. 


Dynamics  and  Control  of  a  Tethered  Spacecraft-  A  Brief  Overview 

Has  some  history  of  the  idea  of  EDTs  and  past  missions  (Gemini).  Description  of  model  and  various 
control  schemes,  including  tension  control,  offset  control,  etc.  See  also  AIAA  1991-1002  below. 


Modi  On  the  Control  of  Tethered  Satellite  System 

Dynamics  of  tether.  Three  different  LQR  controllers  (Thrusters,  tension  and  offset  control)  for 
stationkeeping  and  retrieval.  See  also  AIAA  1990-1 195  prev  entry. 

Mankala  Equilibrium-to-Equilibrium  Maneuvers  of  Rigid  EDTs 

■545  Note  discusses  a  2-D  in  plane  libration  stability  of  a  varying  resistance  EDT.  Tether  modeled  as  rigid 

rod  in  equitorial  circular  orbit.  Feedback  linearization  is  used  to  provide  the  control  history  to  move 
from  one  equilibrium  position  to  another.  Stability  is  not  really  addressed,  but  phase  plot  of  the 
model  used  shows  stable  focus  for  a  given  set  of  tether  params.  Equilibrium  points  are  expressed  in 
terms  of  L  and  radial  distance,  r.  Interesting. 


Mankala  Equilibrium-to-Equilibrium  Maneuvers  of  Flexible  EDT  in  Equatorial  Orbits 

Discusses  shape  of  massless  flexible  EDT.  Control  resistor  on  a  flexible  massless  EDTw/  dominant 
end  mass  in  equatorial  orbit  (B  field  is  perp  to  orbit  and  only  er-et  plane)  Tether  takes  shape  of  an  arc 
of  circle  and  tether  tension  is  constant.  Equilibria  are  unstable  w/o  controller,  but  are  stabilizable 
using  either  linear  controller,  non-linear  controller  or  combo  of  both. 


Williams  Libration  Control  of  Flexible  Tethers  Using  Electromagnetic  Forces  and  Movable  Attachment 

See  also  JGCD  v.27  n5  2004.  Good  ref  list  and  what  is  in  them.  EOM  included-  rigid  and  flexible, 
no  drag,  mother  satellite  mass  dominant,  circular  orbit.. 

Fujn  Nonlinear  Dynamics  of  Tethered  Subsatellite  system  During  Stationkeeping  Phase 
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AIAA  2001-1141 


AIAA  2001-3980 


AIAA  1759-102 


JGCD 


JGCD  1989  0731-5090  Vol 
12#3  (431-433) 


book 


book 


JAS  Vol  48#4  2000  p449-476 


JCGD  v27  n5  2004 


AIAA  3565-487 
Journal  of  Spacecraft  and 
Rockets,  2000,  Vol37#2  187- 
196 


JoVibeandAcousticsl  27_2_20 
05 

Journal  of  Vibration  and 
Acoustics,  Col  127,#2pp144- 
156 

Eur  Jour  of  Mech 

Vol  9,  #2,  1990,  pp207-224 

AIAA-4 147-252 


Poincare  method  used  demonstrating  chaotic  behavior.  Uses  Lyaponov  exponents  and  generates 
bifurcation  maps.  Models  dominant  mass  w/  mt->0. 

Van  Noord  EDT  Tape  Tether  Performance  in  LEO 

Tether  survivability  w/  wider  and  shorter  tethers.  Reduce  Drag,  increase  life.  Includes  twisting  in 
model.  Charts  show  decay  time  vs.  tether  width.  Shows  boosting  by  tether  of  various  widths. 


Van  Noord  EDT  Optimization  for  STEP-AirSEDs  Mission 

Design  optimization  considering  survivability,  drag,  current  collection,  thrust  produced,  tether 
strength,  thickness,  etc.  Has  charts  showing  tether  sever  risk  vs.  width.  Orbit  transfer  time  vs.  width. 


Effect  of  Electromagnetic  Forces  on  Orbital  Dynamics  of  Tethered  Satellites 
JGCD2005-G05-162 


Misra  Breakwell  Memorial  Lecture-  Dynamics  and  Control  of  Tethered  Satellite  Systems  (Sept  29,  2003) 

Only  owned  by  SISTI.  OCLC#57023082.  54th  international  Astro  Congress.  Can't  get  my  hands  on 
this  one. 


Von  Flotow  Some  Approximations  for  the  Dynamics  of  Spacecraft  Tethers 

Explains  why  simple  tether  model  is  good  enough.  Walks  through  methodical  approach  to 
approximate  dynamic  modeling.  Discusses  curvature  and  stress,  strain  relations. 


Misrah  Comments  on  "Some  Approximations  for  the  Dynamics  of  Spacecraft  Tethers" 

Misrah  took  issue  w/  some  of  the  assumptions  in  Von  Flotow's  paper.  A  followon  paper  by  Von 
Flotow  takes  issue  with  Misrah's  issue.  I'll  stay  out  of  it,  but  the  issue  only  relates  to  deployment,  not 
station-keeping  dynamics. 

Effects  of  Atmospheric  Density  Gradient  on  Control  of  Tethered  Subsatellites 
Need  to  obtain  for  longer  tethers 

Beletsky  Dynamics  of  Space  Tether  Systems 

Advances  in  the  Astronautical  Sciences  (An  AAS  publication),  Vol  83.  Checked  out  from  Library. 

Penzo  Tethers  in  Space  Handbook 

NASA  report  edited  by  Paul  A.  Penzo. 

Palaez  A  New  Kind  of  Dynamic  Instability  in  Electrodynamic  Tethers 

ED  tether  modeled  as  a  rigid  rod  w/  point  endmasses.  Constant  tether  current  causes  constant  energy 
being  pumped  into  system  causing  instability  w/  current  on.  There  are  no  equilibrium  positions  (circ 
inclined  orbit).  Equations  have  periodic  solutions.  Does  not  consider  variable  current,  tether 
flexibility  or  damping.  Floquet  theory  used  for  periodic  solutions  and  stability  analysis. 


Williams  Libration  Control  of  Flesxible  Tethers  Using  Electromagnetic  Forces  and  Movable  Attachment 

Good  ref  list  and  what  is  in  them.  EOM  included-  rigid  and  flexible,  no  drag,  mother  satellite  mass 
dominant,  circular  orbit.  See  also  AIAA  2003-5781 

Linear  Stability  Analysis  of  EDTs 
JGCD  G05-109 

Forward  Terminator  Tether:  A  Spacecraft  Deorbit  Device 
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JGCDVol  8,  #1,  Jan-Feb 
1985 


Mankala 


Dynamic  Modeling  and  sumulation  of  Satellite  Tethered  Systems 
Models  tether  shape  dynamics. 


AIAA-2947-955 


Matteis 


Ketchichian 


Wiesel 


Stevens 


Equilibrium  of  a  Tether-Subsatellite  System 

Resonance  due  to  aero  gradient  forces  on  subsatellite  can  cause  instability.  Sensitive  to  atm  model. 
Assumptions-  tether  has  no  mass  or  aero  forces  on  it,  no  bending.  Alt  ~1 10  km.  Assuming  relative 
constant  atmosphere,  no  problem. 

Trajectory  Optimisation  Using  Non-Singular  Elements  and  True  Longitude 
Contains  eom  w/  states  that  avoid  singularities 

Optimal  Many-Revolution  Orbit  Transfer 
Multiple  time-scale  problem. 

Preliminary  Design  of  Earth-Mars  Cyclers  Using  Solar  Sails 
Optimal  control  methods  for  low  thrust  orbital  maneuvering 
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that  is  solved  using  a  pseudospectral  method.  Optimal  control  solutions  are  determined  that  maneuver  an  electrodynamic  tether  to  new  orbits  over  long  time 
scales  while  managing  librational  motion  using  only  current  in  a  wire.  Results  are  simulated  using  a  higher  fidelity  “truth"  model  to  validate  the  controller 
performance. 
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